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Abstract 


The  interaction  of  water  waves  with  an  ice  cover  in  river  channels  is  studied  using 
one-dimensional  governing  equations  which  describe  fluid  continuity,  momentum,  and  the 
ice  cover  response.  The  ice  cover  is  assumed  to  be  a  continuous,  thin  elastic  plate.  Both 
linearized  and  nonlinear  forms  of  the  above  three  equations  are  analyzed  in  this  study. 

The  effects  of  a  constant  compressive  axial  force  along  an  ice  cover  on  the  propagation 
characteristics  of  linear  water  waves,  i.e.  celerity,  attenuation,  and  group  velocity,  are  first 
investigated  over  the  entire  spectrum  of  wavelengths  on  the  basis  of  linear  stability  theory. 
It  is  found  that  only  when  an  axial  force  approaches  a  critical  value  and  the  wavelength  is 
around  2irl ,  where  l  is  the  characteristic  length  of  the  cover,  the  effect  of  the  axial  force  is 
significant. 

The  nonlinear  terms  in  the  water  equation,  which  are  ignored  in  the  linear  theory, 
become  significant  in  highly  unsteady  flows  such  as  surges  caused  by  ice  jam  releases.  An 
analytical  solution  for  the  nonlinear  shallow  wave  equations  is  obtained.  The  effects  of  ice 
inertia  and  an  axial  force  on  the  nonlinear  waves  are  examined. 

The  stresses  induced  in  ice  cover  by  both  linear  and  nonlinear  waves  obtained  from  the 
above  mentioned  analysis  are  used  to  investigate  the  possible  formation  of  closedly-spaced 
transverse  cracks  often  observed  during  ice  cover  breakups.  The  formation  of  closedly- 
spaced  transverse  cracks  is  the  key  mechanism  for  river  ice  breakup  and  the  initiation 
of  breakup  ice  runs.  The  validity  of  the  linear  theory  is  found  to  be  very  limited.  The 
nonlinear  analysis  provides  results  in  close  agreement  with  field  observations.  These  results 
provide  basic  explanation  of  the  formation  mechanism  of  transverse  cracks  which  initiate 
the  breakup  ice  runs. 
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Abstract 


The  interaction  of  water  waves  with  an  ice  cover  in  river  channels  is  studied  using 
one-dimensional  governing  equations  which  describe  fluid  continuity,  momentum,  and  the 
ice  cover  response.  The  ice  cover  is  assumed  to  be  a  continuous,  thin  elastic  plate.  Both 
linearized  and  nonlinear  forms  of  the  above  three  equations  are  analyzed  in  this  study. 

The  effects  of  a  constant  compressive  axial  force  along  an  ice  cover  on  the  propagation 
characteristics  of  linear  water  waves,  i.e.  celerity,  attenuation,  and  group  velocity,  are  first 
investigated  over  the  entire  spectrum  of  wavelengths  on  the  basis  of  linear  stability  theory. 
It  is  found  that  only  when  an  axial  force  approaches  a  critical  value  and  the  wavelength  is 
around  2tt/,  where  l  is  the  characteristic  length  of  the  cover,  the  effect  of  the  axial  force  is 
significant. 

The  nonlinear  terms  in  the  water  equation,  which  are  ignored  in  the  linear  theory, 
become  significant  in  highly  unsteady  flows  such  as  surges  caused  by  ice  jam  releases.  An 
analytical  solution  for  the  nonlinear  shallow  wave  equations  is  obtained.  The  effects  of  ice 
inertia  and  an  axial  force  on  the  nonlinear  waves  are  examined. 

The  stresses  induced  in  ice  cover  by  both  linear  and  nonlinear  waves  obtained  from  the 
above  mentioned  analysis  are  used  to  investigate  the  possible  formation  of  closedly-spaced 
transverse  cracks  often  observed  during  ice  cover  breakups.  The  formation  of  closedly- 
spaced  transverse  cracks  is  the  key  mechanism  for  river  ice  breakup  and  the  initiation 
of  breakup  ice  runs.  The  validity  of  the  linear  theory  is  found  to  be  very  limited.  The 
nonlinear  analysis  provides  results  in  close  agreement  with  field  observations.  These  results 
provide  basic  explanation  of  the  formation  mechanism  of  transverse  cracks  which  initiate 
the  breakup  ice  runs. 
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Chapter  1 


Introduction 

1 . 1  Background 

The  breakup  of  a  river  ice  cover  can  generally  be  classified  as  a  thermal  breakup  or  a 
mechanical  breakup.  A  thermal  breakup  is  essentially  a  thermal  meltout  process  with 
the  floating  ice  cover  deteriorating  and  melting  in  place  with  insignificant  ice  movement. 
This  type  of  breakup  occurs  when  the  river  discharge  remains  relatively  steady  during 
the  spring  breakup  period.  The  physics  of  thermal  breakup  is  relatively  well  understood. 
A  mechanical  breakup  is  the  fragmentation  of  a  floating  river  ice  cover  by  hydraulic  and 
mechanical  forces  associated  with  the  ice  rubble  load  and  the  changes  in  river  discharge 
and  water  level.  The  driving  forces  that  contribute  to  a  mechanical  breakup  include:  1) 
frictional  forces  due  to  water  current  and  wind;  2)  streamwise  component  of  the  weight  of 
the  cover;  3)  the  longitudinal  and  vertical  forces  acting  on  the  ice  cover  by  the  ice  rubble 
in  an  ice  run  from  upstream;  and  4)  pressure  variation  on  the  underside  of  the  ice  cover 
produced  by  the  water  wave. 

A  mechanical  breakup  often  leads  to  severe  ice  runs.  Breakup  ice  runs  and  the  associated 
ice  jams  can  be  destructive  to  hydraulic  structures  and  shoreline  properties.  In  addition  to 
these  concerns,  the  problem  of  mechanical  breakup  is  also  important  for  winter  operations 
of  dams  and  hydropower  stations,  where  the  flow  regulation  is  limited  by  the  stability  of 
the  ice  cover. 

In  the  last  ten  years,  significant  progress  has  been  made  on  the  understanding  of  the 
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dynamics  of  ice  runs  and  ice  jams  (Shen  et  al.  1990,  1993,  Liu  and  Shen  1998).  However, 
since  a  clear  understanding  of  the  mechanics  of  river  ice  breakup  is  not  available,  predictions 
of  the  occurrences  of  breakup  ice  runs  and  ice  jams  still  can  not  be  made.  The  lack  of 
understanding  the  mechanics  of  river  ice  breakup  is  mainly  due  to  the  lack  of  understanding 
on  the  phenomena  of  fracture  of  an  ice  cover  under  the  action  of  surge  waves,  which  is  the 
key  to  the  initiation  of  ice  runs  from  a  floating  ice  cover. 

1.2  Mechanical  Breakup 

Prowse  and  Demuth  (1989)  categorized  mechanical  breakup  into  pre-frontal  and  frontal 
modes.  During  the  initial  stage  of  the  breakup,  pre-frontal  mode  represents  the  breakup  of 
an  ice  cover  into  large  sections  of  ice  sheets.  Frontal  mode  represents  failure  coincident  with 
the  surge  wave  and  ice  runs  during  the  later  stage  of  the  breakup.  Early  in  spring,  shortly 
after  the  beginning  of  spring  runoff,  uplift  pressures  develop  on  the  underside  of  shore- 
fast  ice  cover.  Longitudinal  cracks  form  along  the  river.  Based  on  the  classical  theory 
of  beams  on  an  elastic  foundation,  Billfalk  (1981)  and  Beltaos  (1985,  1990)  developed 
theoretical  explanations  for  the  formation  of  longitudinal  cracks  along  the  river.  These 
analyses  assumed  that  the  dynamic  effects  are  negligible,  and  the  longitudinal  gradient  of 
the  uplift  pressure  is  small. 

The  formation  of  longitudinal  cracks  is  a  prelude  to  the  breakup.  It  is  often  followed 
by  the  formation  of  far-spaced  transverse  cracks.  Shulyakovskii  (1972)  and  Beltaos  (1984) 
showed  that  in  meandering  channels  transverse  cracks  can  form,  spaced  in  the  order  of  1000 
times  the  cover  thickness,  due  to  the  bending  moment  produced  by  the  in-plane  hydraulic 
forces.  Beltaos  (1990)  defined  the  onset  of  breakup  at  a  given  site  as  the  time  when  the 
local  ice  cover  is  set  in  motion.  He  derived  the  criterion  for  the  onset  of  breakup  based  on 
the  boundary  constraints  of  the  river  plane  geometry  in  allowing  the  movement  of  the  large 
ice  sheets  formed  by  far-spaced  transverse  cracks.  However,  these  movements  are  limited 
and  usually  do  not  lead  to  ice  runs  until  closedly-spaced  transverse  cracks  form  when  a 
surge  wave  passes  by. 

The  theories  of  longitudinal  and  transverse  crack  formations  discussed  above  may  be 
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used  to  describe  the  pre-frontal  modes.  The  frontal  modes  occur  during  the  later  stage  of 
the  breakup  period,  and  leads  to  ice  runs,  which  are  caused  by  a  rapid  increase  of  river 
discharge  either  due  to  runoff  increase  or  due  to  an  ice  jam  release  from  upstream.  The 
process  is  much  more  dynamic,  and  may  not  be  analyzed  with  static  formulations  discussed 
above.  Parkinson  (1982)  provided  the  following  vivid  account  on  the  dynamics  of  ice  cover 
fracture  during  the  breakup  event  on  the  Liard-Mackenzie  River  from  a  surge  wave  acting 
on  a  competent  ice  cover: 

Upon  the  arrival  of  the  flood  wave,  there  was  a  rapid  rise  of  the  cover.  ...Cracks 
spaced  at  50  to  200  m  apart  appeared  across  the  ice  cover,  followed  by  small 
ridges  of  crushed  ice  formed  along  the  cracks.  Beyond  this  stage,  ...  the  water 
level  rise  continued  and  the  cover  re-adjusted  and  moved  slowly  from  10  to  100 
m,  until  it  wedged  against  the  shore  again.  This  movement  was  accompanied  by 
widespread  breakage  of  the  cover,  with  crack  lines  and  associated  ridging  spaced 
at  10  to  50  m.  This  in-situ  fracture  of  the  ice  cover  when  the  front  of  the 
surge  wave  was  passing  occurred  over  long  reaches  of  the  river.  The  broken  ice 
front  was  moving  at  a  speed  approximately  equal  to  that  of  the  free  surface  surge 
wave.  The  broken  ice  did  not  move  downstream,  and  the  water  level  returned  to 
approximately  the  pre- flood  level  following  the  passage  of  the  wave —  Prolonged 
ice  run  occurred  only  after  the  water  level  rose  a  couple  days  later  to  the  point 
when  the  broken  ice  could  be  lifted  and  carried  downstream. 

This  description  and  similar  observations  by  Gerard  et  al.  (1984)  and  Prowse  (1986) 
clearly  show  the  importance  of  the  effect  of  a  surge  wave  on  the  formation  of  breakup  ice 
runs.  In  their  attempts  to  explain  the  formation  of  closedly-spaced  transverse  cracks  due 
to  a  surge  wave  similar  to  that  described  above,  Billfalk  (1982)  and  Beltaos  (1985)  assumed 
static  pressure  distribution  on  the  underside  of  the  ice  cover  to  represent  the  uplift  due  to 
surges  and  found  that  a  water  surface  slope  larger  than  0.005  is  required  to  break  the  ice 
cover.  Such  a  large  slope  can  not  be  found  in  rivers  with  an  ice  cover.  Daly  (1995)  used 
linear  analysis  to  study  the  interaction  between  river  waves  and  ice  cover,  and  suggested 
the  possibility  of  the  formation  of  cracks  spaced  at  10  m  or  less  by  waves  in  the  gravity 
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wave  range  with  small  wave  amplitude.  This  phenomenon  was  not  observed  in  the  field. 

In  this  study,  the  linear  analysis  of  Daly  (1993,  1995)  is  extended  to  include  the  effect 
of  an  axial  compressive  force.  The  limitation  of  linear  analysis  is  examined.  A  nonlinear 
analysis  is  carried  out  on  the  interaction  of  a  floating  river  ice  cover  with  shallow  water 
waves,  and  the  formation  of  transverse  cracks  during  the  passage  of  the  waves. 


4 


Chapter  2 


Governing  Equations 

It  is  commonly  accepted  that  the  effect  of  ice  cover  on  water  wave  propagation  in  river 
channels  is  limited  to  the  additional  flow  resistance  due  to  the  presence  of  an  additional 
rough  boundary  and  the  reduction  in  flow  cross  section  due  to  the  ice  cover  submergence. 
Hydrostatic  condition  is  considered  to  be  valid.  However,  for  a  rapidly  varying  channel 
flow,  the  hydrostatic  assumption  may  not  be  valid.  The  wave  may  be  modified  by  the 
inertia  and  resistance  to  bending  of  the  cover. 

The  governing  equations  for  one-dimensional  unsteady  flow  in  a  wide,  rectangular  ice- 
covered  channel  is  presented  in  this  chapter.  The  ice  cover  is  assumed  to  be  an  intact, 
uniform,  thin  elastic  plate.  The  ice-cover  density  and  thickness  are  assumed  to  be  constant 
in  space  and  time.  The  fluid  density  changes  in  response  to  changes  of  fluid  pressure. 


2.1  Water  Mass  Conservation  Equation 

Considering  a  control  volume  in  a  one-dimensional  channel  as  shown  in  Fig.  2.1,  the  mass 
conservation  of  water  in  the  differential  control  volume  is  given  as: 
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Figure  2.1:  Differential  element  for  mass  conservation  equation 


\{pAu)x  -  (pAu)x+&x]At  =  (pAAx)t+At  ~  ( pAAx)t 
In  differential  form,  this  equation  becomes 


djpA )  djpAu) 

dt  dx 


1  Dp  IDA  du 
~p~Dt  +  A  Dt  +  dx 


(2.1) 


where  -§-t  =  total  derivative;  x  =  longitudinal  distance;  t  =  time;  A  =  channel  cross- 
sectional  area;  u  =  cross  section-averaged  flow  velocity;  p  =  fluid  density.  If  the  flow  is 
incompressible,  Eq.  2.1  reduces  to  the  traditional  continuity  equation  for  open  channel 
flows. 

Assuming  that  changes  in  fluid  density  are  caused  only  by  pressure  deviations  P'  from 
hydrostatic,  thus 


IDp  _  1  DP' 
p  Dt  Ew  Dt 


(2.2) 
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where  Ew  =  bulk  modulus  of  elasticity  of  water  with  a  value  of  about  2.04  x  10 9N/m2 
near  0 °C.  Equations  2.1  and  2.2  give 


1  £>P'  1  DA  du 

Ew  Dt  +  A  Dt  dx 


(2.3) 


For  a  wide  uniform  rectangular  channel, 

1  DP'  1  Dd  du  n 

- 1 - 1 - =  0 

Ew  Dt  d  Dt  dx 

where  d  =  depth  of  flow  from  channel  bottom  to  ice  cover. 


(2.4) 


2.2  Water  Momentum  Equation 

Considering  a  control  volume  of  length  Ax  as  shown  in  Fig.  2.2,  the  resultant  force  on  the 
control  volume  in  the  longitudinal  direction  is 


F\  —  F2  +  Wx  —  Ff  —  — '7  A  xA 


d  ( J  Pi  P  \  T 

—  I  d  +  z  4 - 77  4  14  — 

dx  \  p  pg)  7  R 


(2.5) 


in  which  R  =  j,  hydraulic  radius;  r  =  average  shear  stress  due  to  the  bed  and  ice  cover; 
rj  =  ice  cover  thickness;  7  =  unit  weight  of  water.  The  term  within  the  parenthesis  is  the 
piezometric  head: 


H  =  z  +  d+  —r)  +  —  (2.6) 

P  P9 

The  right-hand  side  terms  of  Eq.  2.5  can  be  derived  as  the  following.  The  pressure 
forces  acting  on  the  control  volume  are 


Fi  =  APi  +  jAz  =  AP[  4-  gPi%A  +  jAz 


F2  =  AP2  +  jAz  4-  7  A  7^  =  AP'2  4-  g  Pity  A  +  7A2  4-  7 


where  z  =  centroid  of  cross  section  A.  The  weight  of  water  in  the  longitudinal  direction  is 
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Figure  2.2:  Definition  sketch  for  momentum  equation 


Wx  =  jAAx  sin  0  —  —jAAx 


dz 

dx 


where  0  «  1.  The  total  friction 


Ff  =  ( Tbpb  +  Tipi}  Ax  =  rpAx 


where  p  =  wetted  perimeter.  According  to  Reynolds  transport  theorem: 


Q  rx+Ax 

Fx  =  diJx 
\9{pA)  d{pAa) 

[  dt  dx 

.  .  ,du  du 

= pA^ai + “&> 


puAdx  +  (pAu2)  2  —  {pAu2)  i 


du  du 

uAx  +  pAAx^+u-^) 


Defining  the  friction  slope  as 


„  (npt  +  npi ) 

S,=  7A 


The  momentum  equation  becomes 


du 

dt 


du  dH  _  n 
+sS/=0 


(2.8) 


When  the  vertical  acceleration  of  the  flow  is  significant,  and  the  flow  is  incompressible, 
Eq.  2.8  becomes  (Chaudhry  1993) 


du  du  dH  _  d? ,  d3u  d3u 

»+"&+9&+9S/  =  T  + 


(2.9) 


The  Boussinesq’s  term  on  the  right-hand  side  of  Eq.  2.9  represents  the  effect  of  vertical 
acceleration,  when  the  flow  is  not  hydrostatic. 

For  a  rectangular  channel,  the  friction  slope  becomes 

2f 


Sf  = 


7  d 
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in  which  f  =  |(r,  +  rfc)  =  |pu2,  and  /  =  average  friction  factor  of  the  bed  and  the  ice 
cover.  Noting  that  r]  and  z  do  not  vary  with  time,  then  the  mass  conservation  equation 
(2.4)  becomes 


(2.10) 


where  a  =  Since  the  celerity  of  acoustic  wave  cQ  = 


~  anc* 


celerity  of  gravity  wave  c  =  y/gd,  the  parameter  a  can  be  approximated  by  a  ~  %. 

If  a  ~  0,  i.e.  c  -C  cQ  ,  the  incompressible,  hydrostatic  flow  condition  prevails,  and  Eq. 
2.10  reduces  to  the  conventional  momentum  equation  for  open  channel  flow  with  a  floating 
cover. 


Dd  du 
Dt  dx 


=  0 


(2.11) 


2.3  Ice  Cover  Response  Equation 

The  ice  cover  is  assumed  to  be  a  homogeneous  thin,  continuous  elastic  plate.  A  differential 
element  of  ice  cover  per  unit  width,  as  shown  in  Fig.  2.3,  is  subjected  to  the  hydrodynamic 
pressure,  P,  on  the  underside  of  the  cover  from  the  underlying  water,  in  addition  to  water 
drag,  wind  drag  and  gravity.  The  cover  is  also  subjected  to  an  axial  force  resulting  from  the 
combination  of  fluid  drag  force,  the  component  of  the  gravitational  force  in  the  longitudinal 
direction  resulting  from  the  slope  of  the  water  surface,  and  the  loading  at  the  upstream  end 
of  the  intact  cover  due  to  ice  rubble  and  fractured  ice  cover.  This  axial  force  is  partially 
balanced  by  the  intermittent  bank  resistance  on  the  cover  along  the  longitudinal  shore 
cracks. 

The  momentum  equation  for  the  ice  element  in  the  x-direction,  neglecting  longitudinal 
inertia,  is 

OF 

F(x)  -  (F(x)  +  As)  +  A xtx  +  Axrjpig  sin  9  =  0 
ox 

or 

dF  . 

—  =Tx  +  rjpigsme 
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Figure  2.2:  Definition  sketch  for  momentum  equation 


where  rx  is  the  x-component  of  r  =  tw  +  Ta  —  -p,;  tw  and  ra  are  shear  stresses  on  ice  due 
to  current  and  wind,  respectively;  and  T(,  =  tb/B  with  tb  =  bank  friction  force  per  unit 
cover  length;  and  B  =  cover  width;  F(x)  is  the  compressive  axial  force  per  unit  width.  The 
axial  force  F(x )  may  be  assumed  to  be  a  constant  N  for  a  uniform  channel  with  constant 
ra  and  tw,  except  near  the  leading  edge  of  the  ice  cover. 

The  momentum  equation  for  the  element  in  vertical  direction  is 

V(x)  -  (V{x)  +  — Arr)  ~  Axrjpig  cos  6  +  A xP(x)  =  Ax- 

in  which,  P(x)  is  the  uplift  pressure  at  the  interface. 

For  6  «  1,  this  equation  becomes 

dV  v  d2d  (0 19x 

-fa  +  VPiQ  -  P(x)  =  -PWqI 2  (212' 

in  which,  P(x)-P'{x)  +  gp^g  is  the  net  upward  pressure  on  the  underside  of  the  ice  cover 
at  its  deformed  state.  The  pressure  P(x)  is  a  combination  of  the  hydrodynamic  pressure 
P(x)  and  the  net  buoyancy  (Beltaos  1990). 

Taking  moment  with  respect  to  O  gives: 

-M{x)  +  (M{x)  +  ^Az)  +  (V(*)  +  ^A x)Ax 

-F(x) Ax^  +  gpiijAx^-  -  P(x)Ax ^  =  0 


i.e. 


or 


£  +  v-r(„!-o 


d2M  dV  ,  .d2d  dJFdd_0 

dx2  +  dx  ^dx2  dx  dx 


(2.13) 


For  a  thin  plate,  the  bending  moment  is  related  to  the  curvature  of  the  deflection  curve 
of  the  plate  (Ugural  1995),  i.e.  M  =  For  small  deformations  (|  |<  1),  the 

curvature  ±  =  [1+ff/g ^  =  0  +  (-§)(g)20  +  -  »  0-  Thus  M  can  be  approximated 
by  M  =  —  yt^20-  Assuming  F(x)  as  a  constant  N.  Equations.  2.12  and  2.13  give  the 
following  ice  cover  equation  of  motion: 

£)2  j 

(2.14) 


El  d*d  d2d  ,  d2d  n 

+  PiV^-T(x)  +  N^  =  0 


1  —  v2  dx4 


dx 2 
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or 


M&d  A&d  F(x)  Nd2d 
pg  dt2  3.x4  pg  pg  dx2 

in  which, 


El 

1/4 

i fE 

pg{\  -  v 2) 

12pg(l  -  v2) 

(2.15) 


(2.16) 


The  characteristic  length  of  the  ice  cover,  Z,  represents  the  stiffness  of  the  cover.  This 
equation  is  the  same  as  the  equation  for  elastic  plate  on  an  elastic  foundation  with  modulus 
k  —  pg.  The  elastic  modulus  of  ice,  E,  is  in  the  range  of  0.4  ~  9.8  GPa,  and  the  Poisson’s 
ratio  v  =  0.35  (Daly  1993). 

Equations  2.14  or  2.15  has  been  linearized  based  on  the  assumption 

(g)2  «  1  (2.17) 

This  implies  that  the  nonlinear  term  (f^)2f^r  is  of  an  order  smaller  than  the  term  in 
the  expression  for  the  ice  cover  curvature  £.  The  validity  of  this  approximation  will  be 
discussed  in  Sections  3.2.3,  4.2.2  and  4.2.4. 
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Chapter  3 


Linearized  Model 


In  this  chapter,  wave  propagation  in  ice-covered  channels  will  be  analyzed  using  the  method 
of  linear  perturbation,  in  which  the  waves  are  represented  by  sinusoidal  perturbations.  The 
length  of  waves  can  be  varied  from  very  long  kinematic  flood  waves  to  very  short  pressure 
waves.  The  general  approach  follows  that  of  Ponce  and  Simons  (1997),  who  analyzed  the 
wave  propagation  in  open  channels.  Daly  (1993,  1995)  studied  the  fracture  of  the  ice  cover 
by  river  waves,  using  the  same  approach.  In  Daly’s  analysis,  the  axial  force  acting  along 
the  cover  was  neglected.  The  purpose  of  the  present  chapter  is  to  introduce  the  effect  of 
the  axial  force  acting  along  the  cover  and  to  examine  the  validity  of  the  linear  theory  in 
the  analysis  of  the  breakup  of  ice  cover  due  to  wave-ice  cover  interaction. 

3.1  Linear  Perturbation  Analysis 

The  governing  equations  2.8,  2.10,  and  2.15  must  satisfy  the  unperturbed  steady  uniform 
flow  condition  for  which  u  =  uo,  d  =  do,  H  =  Ho,  and  r  =  To  as  well  as  the  perturbed  flow 
for  which  u  =  u0  +  u,  d  =  d0  +  d\  H  =  H0  +  H' ,  and  r  =  r0  +  r',  where  superscripted 
variables  represent  a  small  perturbation  to  the  steady  uniform  flow.  The  gradients  of  the 
steady  uniform  variables  are  zero  in  x  and  t  except  that 


Since  the  perturbations  in  ice  thickness  r}  and  bed  elevation  z  are  not  allowed,  Ho  — 
z  +  do  +  for  the  hydrostatic  condition.  Equation  2.6  gives  H'  =  H  —  H0  as: 

H'  =  d  +  —  (3.2) 

P9 

Substituting  the  perturbed  variables  into  Eqs.  2.6,  2.9  and  2.12  and  subtracting  the 
undisturbed  equations  yield  the  following  linearized  perturbation  equations  after  neglecting 
higher-order  terms: 


„  x/dd  d(i ,  ,dH' 


atf.  ,  du  „ 

+u°-fc)+d°~te  -° 


du  du'  OH'  ,2u  d 

—  +  u0—  +  g —  +  gS0{—  —  ~r)  =  0 
ot  ox  ox  Uo  do 


piT)  d2d  4  d*d  ,  y.  N  d2d 
-^-dW  +  lW~{H  -d)  +  7gW 


In  order  to  provide  a  convenient  way  to  examine  the  various  wave  models,  Eq.  3.4  is 
recasted  as 


lcdu  auodu  8H’  Ift/2u  n  fo 

JL  +—1  +p  +kS0( - T-)  =  0  (3‘6) 

g  Ot  g  Ox  Ox  u0  d0 

in  which  lc,  a,p ,  k  are  integers  that  can  take  a  value  of  either  0  or  1  ,  depending  on  which 
terms  are  used  to  describe  the  wave  motion. 

If  the  flow  is  assumed  to  be  incompressible,  the  effect  of  including  the  vertical  inertia 
of  water  may  be  analyzed  by  using  a  linearized  Boussinesq’s  approximation.  The  above 
equation  becomes  (Chaudhry  1993): 


lcdu  auQdu  dH '  .2u  d  -  d%  d3u  d3u\ 

g  dt  g  dx  dx  uq  do  3 g  dx2dt  dx 3 
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in  which  b  is  an  integer  that  can  take  a  value  of  either  0  or  1. 

The  propagation  of  shallow  water  waves  is  controlled  by  the  balance  of  the  various  forces 
included  in  the  momentum  equation.  In  Eq.  3.7,  the  first  term  represents  the  local  inertia 
term,  the  second  term  represents  the  convective  inertia  term,  the  third  term  represents  the 
pressure  differential  term,  the  fourth  term  accounts  for  the  friction  and  bed  slopes,  and  the 
fifth  term  represents  the  vertical  inertia  of  water.  Various  wave  models  can  be  constructed, 
depending  on  which  of  these  terms  is  assumed  to  be  negligible  when  compared  with  the 
remaining  terms.  The  wave  models  are:  (1)  Kinematic  wave,  lc  =  a  —  p  =  b  =  0,k  =  1; 
(2)  diffusion  wave,  lc  =  a  =  b  =  0,p  =  k  =  1;  (3)  steady  dynamic  wave,  lc  =  b  =  0,a  =  p  = 

k  =  1;  (4)  gravity  wave,  lc  —  a  —  p  =  1,  k  =  b  —  0;  (5)  dynamic  wave,  lc  —  a  =  p  =  k  =  l, 

6  =  0;  and  (6)  Boussinesq’s  approximation,  lc  =  a  —  p  =  k  =  b  =  \. 

Assuming  the  perturbations  are  in  the  following  exponential  form  : 

—  =  de[i{ai~0t)]  (3.8) 

d0 

—  =  ue^ai~pt)]  (3.9) 

u0 

—  =  (3.10) 

d0 

in  which  d,  u,  and  H  =  dimensionless  amplitudes;  x,  t  =  dimensionless  space  and  time 
coordinates;  a  =  a  dimensionless  wave  number;  and  /?  =  a  dimensionless  complex  wave 
propagation  factor.  These  dimensionless  variables  are  defined  as  the  following: 


(3.11) 

X 

Z=  — 

Ml 

(3.12) 

2  tu0 
~To 

(3.13) 

where  L0  =  a  longitudinal  length  scale;  and  /3  =  Pr+i  fa.  The  imaginary  part  describes 
the  dimensionless  wave  amplitude  de^. 
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By  defining  the  real  part  of  P  as 


£ 


Pr  — 


2ir  Lq 
T  «o 


(3.14) 


the  dimensionless  wave  celerity  c  can  be  written  as 


<7 


(3.15) 


where  c  =  and  the  wave  celerity  c  =  ^,  in  which  L  —  wave  length  and  T  —  wave  period. 

Substituting  (3.7),  (3.8),  and  (3.9)  into  (3.3),  (3.5),  and  (3.6)  results  in  the  following 
set  of  equations: 


iap  2k  +  iF2[(ao  —  lc(3)  +  ba2(cr  —  /?)] 
a(a  —  P)  a 

-1  0 


—k 

(1  -  ot)(a  -  P) 
Y 


H 

u 

d 

=  0 


(3.16) 


where 


and 


Y  =  -7 1oF?A2P2  +  l40a4  +  1  -  N0a2 
L0 


No 

Fr 

Vo 

lo 

b 


N 

pgLa 

Up 

\fgd o 

PiV_ 

P  do 

J_ 

Lo 

Wb=\s® 


(3.17) 


(3.18) 

(3.19) 

(3.20) 

(3.21) 

(3.22) 


Equation  3.16  constitute  a  homogeneous  system  of  linear  equations  in  the  unknowns  H, 
u  and  d.  For  the  solution  to  be  nontrivial,  the  determinant  of  the  coefficient  matrix  must 
vanish.  The  determinant  is  a  fourth-degree  polynomial  of  p  with  complex  coefficients.  The 
wave  dispersion  relation  can  be  obtained  from  the  determinant  matrix: 
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(3.23) 


J2(Pj  +  i(li)Pj  1  =  0 

j= i 

where 

px  =  k(— 2lgcr5a  —  3<r)  +  k(2aN0a 3) 

p2  =  k(  2  +  2al^aA)  4-  k(—2N0cr2a) 

p3  =  k(2XF2aa) 

Pi  =  k(—2XF2a) 

p5  =  0  (3.24) 

and , 

ft  =  pa2  4-  p/qCt6  —  a(F2ll<76a)  —  aF2a2  —  N0aA  (p  —  aaF2)  —  boAF2(all<7A  —  aN0a2  4  1) 

ft  =  (Je  +  a) (F2a)  +  (/c  +  a) (F2^a5a)  +  (lc  +  a)(-N0F2a3a)  4  2 ba3F2(al40a4  -  aN0a2  4  1) 

&  =  aXF4a2a  -  p\F2a2  -  lcF2  -  lc(alAaAF2)  4  lc{aF2N0a2) 

—ba2F2[alloA  —  aN0a 2  +  1  +  a2(— ar]0F2  (-p-)2)] 

Lq 

ft  =  - (o  +  /c) XFAaa  -  2bar]0a3F2(Y-)2 

L o 

g5  =  lcXF*a  +  bar]0a2F2(Y-)2  (3.25) 

4o 

In  Eq.3.23,  /?  can  be  analytically  solved  as  a  function  of  cr,  which  can  be  used  to  analyze 
wave  celerity  and  attenuation  under  ice-covered  conditions. 

A  computer  code  is  developed  to  calculate  the  four  complex  roots  of  Eq.  3.23  using  the 
IMSL  routine  DZPOCC.  Before  doing  this,  Eq.  3.23  is  recasted  by  introducing 

/T  =  £[1  +  («F2)-5]-i  (3.26) 

a 

and  recognizing  that  the  wave  number  a  can  be  defined  by 

^=<r/0  (3-27) 
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In  general,  four  roots  exist  in  Eq  3.23.  Two  of  them  describe  the  propagation  of  free 
waves  in  an  infinite  continuously  supported  plate,  and  are  not  relevant  to  the  present  inves¬ 
tigation.  The  remaining  two  roots  describe  two  waves  propagating  along  two  characteristic 
paths. 

3.2  Analysis  of  Results 

3.2.1  Wave  Celerity 

The  propagation  characteristics  of  various  types  of  water  waves  in  open  channels  were 
discussed  in  detail  by  Ponce  and  Simons  (1977).  The  solution  of  Ponce  and  Simons  (1977) 
can  be  considered  as  a  special  case  of  ice-covered  channels  with  r)  =  l  =  N  =  b  =  0.  The 
propagation  celerity  for  open-water  conditions  as  a  function  of  Froude  number  and  wave 
number  are  shown  in  Fig.  3.1,  which  is  the  same  as  Fig.  1  of  Ponce  and  Simon  (1977). 


Figure  3.1:  Dimensionless  wave  celerity  as  function  of  dimensionless  wave  number  in  open 
water  channel 
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This  figure  shows  that  there  are  three  bands  in  the  wave  number  spectrum  for  the 
dynamic  model:  (1)  a  gravity  band  corresponding  to  large  wave  numbers,  in  which  c  is 
1  +  -r,  independent  of  cr;  (2)  a  kinematic  band  corresponding  to  small  wave  numbers,  in 
which  c  is  a  constant  equal  to  3/2.  and  (3)  a  dynamic  band  corresponding  to  intermediate 
values  of  wave  number,  in  which  c  varies  with  both  o  and  Fr. 


Figure  3.2:  Effect  of  axial  force  on  the  dimensionless  wave  celerity  in  ice-covered  channel. 

Under  ice-covered  condition,  the  variation  of  wave  celerity  is  shown  in  Fig.  3.2.  The 
dimensionless  wave  celerity  is  calculated  by  using  the  first  root  of  Eq.  3.23,  which  represents 
the  ’’primary  wave”  propagating  in  the  downstream  direction  (Ponce  and  Simoms  1977). 

As  pointed  out  by  Daly  (1993),  there  are  five  well-defined  bands  in  the  wave  number 
spectrum  for  ice-covered  channels.  As  shown  in  Fig.  3.2,  the  first  band  from  the  left, 
corresponding  to  small  values  of  dimensionless  wave  number,  was  identified  as  the  kinematic 
band  (K),  which  has  a  wave  celerity  c  =  |,  or  c  =  |«o>  in  which  the  gravity  and  frictional 
forces  dominate.  The  next  band  to  the  right,  corresponding  to  larger  dimensionless  wave 
numbers,  is  the  the  dynamic  band  ( D )  over  which  c  changes  rapidly  and  is  a  function  of 
both  <j  and  Fr,  in  which  the  gravity,  frictional,  and  inertial  forces  dominate.  With  a  further 
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increase  in  the  dimensionless  wave  number,  the  gravity  band  ( G ),  in  which  the  gravity  and 
inertial  forces  dominate,  is  found.  The  wave  celerity  is  only  a  function  of  Fr.  For  the  wave 
propagating  downstream,  c  =  1  4-  j-r,  or  c  =  uq  +  (gdo)1^2.  This  is  analogous  to  the  open- 
water  gravity  wave  celerity.  These  three  bands,  i.e.  kinematic,  dynamic,  and  gravity  bands, 
define  the  quasi-open-channel  range.  The  wave  celerity  is  not  affected  by  the  existence  of 
the  ice  cover,  and  independent  of  variations  in  the  ice  properties.  The  governing  equations 
could  be  simplified  by  dropping  the  ice  cover  equation  and  let  P  =0.  The  ice  cover  can 
be  considered  as  floating  at  hydrostatic  equilibrium  over  the  range  of  quasi-open-channel. 

Immediately  above  the  quasi-open-channel  range,  the  ice-covered  channel  wave  celerity 
begin  to  deviate  from  the  corresponding  open-water  wave  celerity  in  the  gravity  wave 
band,  in  which  the  gravity,  inertial,  and  ice-cover  bending  and  inertial  forces  dominate. 
This  band  can  be  used  to  define  the  limit  of  quasi-open-channel  range.  Daly  (1993)  found 
in  this  band,  which  is  named  as  the  ice-coupled  band  (I),  the  dimensionless  wave  celerity 
deviates  sharply  from  the  open  water  condition  and  is  a  function  of  nondimensional  wave 
number,  the  Froude  number,  and  the  characteristic  length  of  ice. 

A  further  increase  in  the  dimensionless  wave  number  leads  to  the  acoustic  band.  It 
is  only  in  this  band  that  water  is  considered  as  compressible  and  the  wave  celerity  are 
equal  to  the  acoustic  wave  velocity  in  water  as  if  waves  were  propagating  in  a  closed  elastic 
conduit.  The  dimensionless  wave  celerity  in  this  band  is  independent  of  the  dimensionless 
wave  number  and  the  ice  properties  and  is  a  function  of  Fr  and  a  only,  c  =  1  ±  FjJ!/2  ■  The 
wave  celerity  is  c  =  u0  ±  (i^.)1/2  —  u0  +  ca. 

The  propagation  of  wave  in  a  floating  ice  plate  is  governed  by  two  restoring  forces 
produced  by  elastic  bending  of  the  plate  and  the  tendency  of  gravity  to  make  the  upper 
surface  of  the  supporting  water  horizontal.  The  compressive  axial  force  acts  against  these 
restoring  forces.  Figure  3.2  shows  the  effect  of  the  axial  force  for  three  different  values  of 
l0  corresponding  to  E  =  1  GPa,  5.5  GPa,  and  10  GPa.  This  shows  that  an  axial  force  in 
the  ice  plate  may  sharply  decrease  the  dimensionless  wave  celerity  in  the  gravity  band  and 
the  ice-coupled  band.  Addition  figures  for  different  values  of  Fr  are  given  in  Appendix  C 
as  Figs.  C.l  and  C.2.  These  figures  show  that  the  increase  in  the  elastic  modulus  of  the 
ice  cover  can  reduce  the  effect  of  axial  force  on  the  wave  celerity. 
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Dimensionless  Wave  Number  2n  l/L=l0c 


Figure  3.4:  Variation  of  the  wave  celerity  with  axial  forces,  c,  relative  to  that  without  axial 
force,  Co,  at  ^  =  1. 


The  influence  of  the  magnitude  of  the  axial  force  on  wave  celerity  can  be  seen  in  Fig. 
3.3,  as  well  as  Fig.  C.3  to  C.ll  for  different  Fr  and  E  values.  Figure  3.4  shows  that 
this  effect  on  wave  celerity  becomes  significant  only  when  the  magnitude  of  the  axial  force 
reaches  the  critical  axial  force,  iVcr,  defined  in  section  3.2.3.  The  ice  thickness  does  not 
have  considerable  effect  on  the  wave  celerity,  as  shown  in  Fig.  3.5. 


Figure  3.5:  Effect  of  ice  cover  thickness  on  dimensionless  wave  celerity  in  ice-covered  channel 

3.2.2  Wave  Attenuation 

The  wave  attenuation  is  described  by  the  dimensionless  amplitude  de^'1.  The  wave  decay 
can  be  quantified  by  the  logarithmic  decrement  of  the  wave  amplitude  over  a  duration 
It  is  defined  as  5  -  /n(ax)  -  Zn(a0),  where  ax  and  o0  =  amplitudes  at  a  time  period 
apart.  Based  on  Eq.  3.13,  this  leads  to 

S  =  In  =  In  =  £cr  (3.28) 

Figure  3.6  as  well  as  Figs.  C.12  and  C.13  depict  that  the  natural  decrement,  5,  of  wave 
amplitude  is  always  negative  for  the  case  considered  here,  indicating  that  waves  propagating 
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in  ice-covered  channels  are  always  attenuated  as  they  progress.  Over  the  quasi-open-channel 
range,  the  ice-covered  wave  attenuation  is  analogous  to  the  open-water  wave  attenuation. 
The  decrement  6  has  its  minimum  value  at  the  lowest  nondimensional  wave  numbers  in  the 
kinematic  band  and  reaches  a  maximum  value  in  the  dynamic  band.  The  value  of  5  closes 
to  a  constant  —  in  the  gravity  band  and  is  much  larger  than  that  in  kinematic  band. 

Again,  similar  to  the  wave  celerity,  the  wave  attenuation  deviates  sharply  from  the  quasi¬ 
open-channel  attenuation  in  the  ice  influenced  bands.  The  attenuation  of  ice-coupled  waves 
decreases  with  the  increasing  wave  number  and  reaches  a  minimum  at  a  nondimensional 
wave  number  slightly  less  than  the  limit  of  the  acoustic  band.  In  the  acoustic  band,  S  is  a 
constant  that  is  approximately  —  It  is  noted  that  dynamic  and  gravity  waves  will  be 
subjected  to  strong  attenuation,  whereas  kinematic  waves,  ice-coupled  waves  at  higher  wave 
numbers,  and  acoustic  waves  will  be  subjected  to  small  attenuation.  Figure  3.6  shows  that 
the  axial  force  increases  the  wave  attenuation  in  the  transition  zone  between  the  gravity 
band  and  ice-coupled  band. 


Figure  3.6:  Effect  of  axial  forces  on  natural  decrement 
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3.2.3  Ice  Cover  Response 

Assuming  that  the  floating  ice  cover  is  subjected  to  hydrostatic  pressure  only,  the  following 
elastic  plate  equation  can  be  obtained  from  Eq.  2.8  and  2.15. 


Pirid2d  4d*d  j  N  d2d 
pg  dt2  dx*  pg  8x2 


=  0 


(3.29) 


This  is  the  equation  of  elastic  plate  supported  by  a  massless  Winkler  foundation.  The  free 
or  homogeneous  wave  celerity  ,  Ch,  in  the  ice  cover  can  be  obtained  by  substitution  of  Eq. 
3.8  into  Eq.  3.29  as 


This  shows  the  compressive  axial  force  can  reduce  the  homogeneous  wave  celerity.  The 
wavelength  Lmin,  at  which  the  minimum  homogeneous  wave  celerity  occurs,  is  found  from 


T-'min  —  2  Ttl 


(3.31) 


and  the  minimum  homogeneous  wave  celerity  is 


(3.32) 


When  Ch  =  0,  the  minimum  value  of  N  obtained  from  Eq.  3.30  using  the  condition 
=  0  is  again  L  =  2ixl.  This  minimum  value  of  N  is 

"i  i  J 


\cH=0  =  Zpgl2  =  2 \fpgEI  =  Na 


(3.33) 


which  is  the  critical  buckling  load  of  the  infinite  beam  (Hetenyi  1955).  Equation  3.32 
shows  that  Cumin  decreases  with  increasing  N  and  Cumin  becomes  zero  when  N  =  Ncr. 
i.e.  N0  =  2 ll 

It  is  well  known  that  when  N  >  Ncr  the  beam  shape  is  not  straight  i.e.  the  assumption 
used  for  Eq.  2.14  is  not  satisfied.  Hence  Eqs.  2.15  and  3.29  are  not  applicable  (Kerr  1972). 
For  this  case  the  necessary  analysis  is  very  involved  and  is  beyond  the  scope  of  this  study. 
In  the  following  analysis,  we  will  only  consider  conditions  with  N  <  Ncr. 
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Figure  3.7:  Dimensionless  free-wave  celerity  C;j  and  dimensionless  wave  celerity  c  in  ice- 
covered  channel  as  functions  of  dimensionless  wave  number 
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The  ice  cover  response  is  described  by  Eq.  3.5.  By  substituting  Eqs.  3.8  and  3.10  for 
the  dynamic  pressure  P'  —  pg(H'  —d)  into  Eq.  3.5,  the  relationship  between  the  amplitude 
of  the  propagating  wave  H  and  the  response  of  the  ice  cover  d  can  be  obtained: 


(3.34) 


(3.35) 


Since  +  (^j)2]  ^  2,  and  reaches  its  minimum  value  of  2  when  L  =  Lmin  =  2nl, 

the  maximum  ice  response  d  is  approached  at  the  wavelength  of  L  =  2irl. 

Equation  3.32  can  be  rewritten  as 


■ + (&)4  - 1  (¥)2]  [i  - 


(3.36) 


This  shows  that  when  c  =  CH,  the  ice  cover  response  d  approaches  infinite.  The 
assumption  of  Eq.  2.14  and  the  small  wave  amplitude  assumption  will  be  violated.  When 
the  axial  force  is  zero,  Fig.  3.7  as  well  as  C.14  to  C.18  showed  that  c  <  C//  over  the  entire 
spectrum  of  possible  wavelength,  although  c  closely  approaches  Ch  in  the  ice-coupled  band 
at  wave  numbers  just  below  the  acoustic  band.  Because  of  the  resolution  of  Fig  3.7,  it  may 
appear  that  c  and  Ch  are  equal  in  this  range,  as  pointed  out  by  Daly  (1993).  The  cover 
response  d  approaches  H  when  the  wave  length  L  is  much  larger  than  2ixl,  and  approaches 
zero  when  L  is  much  less  than  2irl.  When  the  axial  force  is  not  zero  and  N  <  Ncr,  c  is 
still  less  than  CH -  When  N  -»  N^,  c  and  CH  approach  zero  at  wavelength  L  =  2i d  as 
shown  in  Fig.  3.7,  and  d,  approaches  infinity.  This  shows  the  possibility  of  cover  fracture 
at  a  wave  length  2? rl  as  N  approaches  Ncr.  However,  the  linear  perturbation  analysis  may 
not  be  valid  as  d  approaches  infinity. 

To  further  the  analysis,  the  bending  stress  in  the  ice  cover  is  examined.  The  bend¬ 
ing  stress  Sx  in  the  ice  cover  at  a  distance  z  from  the  neutral  axial  produced  by  waves 
propagating  in  the  longitudinal  direction  is 
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(3.37) 


0  Ez  d2d 

Sx~  1  -  1/2  dx2 

The  maximum  bending  stress  in  the  ice  cover  SXtTnax  is  at  the  top  and  the  bottom  of 
the  ice  sheet,  z  —  ±§. 


E  r)d2d 
1  —  v1 2  dx2 


(3.38) 


By  substituting  Eq.  3.8  for  d  and  then  relating  d  to  H  by  Eq.  3.35,  the  maximum 
stress 


HdpEfi  L/2ir\2  (L\2 

X’max  1  —  v2  2  \L)  \2ttJ 


2  N_ 
P9 


i-lr 


-i-l 


1- 


cl , 


(3.39) 


Equation  3.39  and  3.30  gives 


max  — 


Hd0E  t;  1 

[«•(*)’+(£)*- 5 


V* 


mr  2 

99  ° 


(3.40) 


For  an  ice  cover  with  a  maximum  flexural  strength  Sxb  >  the  minimum  wave  amplitude  to 
cause  cover  crack  at  any  wavelength  is  then  obtained  as 


^min 


^0  Hmin{L) 


25x6(1  -I/2)  ,2 
Ef] 


(3.41) 


where  ^  The  amplitude  amin(L)  is  shown  in  Fig.  3.8  as  well  as  C.19  to  C.23 

as  a  function  of  the  dimensionless  wave  number  =  l0<y.  These  figures  show  that  when 
the  axial  force  increases,  the  amplitude  of  the  propagating  wave  required  to  fracture  the 
cover  decreases  when  the  wave  length  is  in  the  vicinity  of  2irl.  The  amplitude  amjn(L)  has 
a  well-defined  global  minimum  value  at  the  wavelength  Lmin  —  2x1.  In  addition,  amtn  (L) 
rapidly  approaches  zero  when  the  axial  force  approaches  the  critical  buckling  load  Ncr  for  a 
homogeneous  wave.  This  means  that,  under  the  assumptions  of  linear  perturbation  theory, 
when  the  ice  cover  is  subjected  to  the  critical  axial  force,  it  will  be  fractured. 
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Considering  that  the  axial  force  is  related  to  the  length  of  ice  cover,  D ,  over  which  there 
is  no  bank  support,  then  Eq.  3.33  gives 

Ncr  =  2  pgl2  =  i DdopgSo 

Using  a  value  of  /  =  5  m  from  Fig.  4.6,  So  =  5  x  10~4,  do  =  5  m,  then  D  =  40  km.  This 
shows  that  the  axial  force  rarely  reaches  the  critical  value  in  the  field  condition  for  natural 
river  ice  cover.  Without  an  axial  force,  the  wave  amplitude  required  to  cause  an  ice  cover 
to  crack  is  O(O.lm)  at  the  wavelength  of  2iri  (Daly  1995).  Hence,  for  both  cases  with  or 
without  axial  load  the  term  (|^)2  1,  and  the  small  deformation  assumption  used  in  Eq. 

2.14  is  satisfied. 


3.2.4  Effect  of  Boussinesq’s  term 


The  one-dimensional  depth-averaged  momentum,  Eq.  2.8  or  3.4  are  based  on  the  assump¬ 
tion  of  hydrostatic  pressure  distribution.  When  the  wave  number  is  large,  the  vertical 
acceleration  of  the  flow  may  be  significant.  Since  the  preceding  discussions  showed  that 
the  fracture  of  an  intact  ice  cover  will  occur  at  a  wave  length  of  27tZ,  the  compressibility  of 
water  associated  with  the  acoustic  band  need  not  be  considered.  This  assumption  is  further 
supported  by  the  fact  that  longitudinal  cracks,  formed  long  before  the  fracture  of  cover  by 
propagating  waves,  enable  the  release  of  excess  pressure  under  the  cover.  The  effect  of  the 
vertical  inertia  of  fluid  may  be  expressed  by  a  linearized  Boussinesq  approximation  as  in 
Eq.  2.9  or  3.7.  Steffler  and  Hicks  (1994)  analyzed  linearized  Boussinesq’s  equation  and 
discussed  the  effect  of  vertical  water  inertia  on  wave  celerity  when  there  is  no  axial  force. 
When  lc  =  a  —  p  =  k  =  b  —  1  is  chosen  in  Eq.  3.7,  the  effect  of  the  Boussinesq’s  term  is 


considered,  i.e. 

1  du  u0du  dH'  ,2u  d  d\  d3u  d3u 

g  dt  g  dx  dx  «o  do  3  <7  ox2ot  ox3 

As  shown  in  Figs.  3.9,  C.24  and  C.25  the  vertical  inertia  of  the  water  causes  a  reduction 
in  celerity.  However,  it  does  not  reduce  minimum  wave  amplitude  required  to  fail  the  ice 


cover,  as  shown  in  Fig.  3.10. 
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Figure  3.9:  (a)Effect  of  Boussinesq’s  term:  without  axial  force,  Steffler  and  Hicks  (1994) 
(b)  with  an  axial  force. 
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Figure  3.10:  Effect  of  vertical  inertia  of  the  water  on  the  minimum  wave  amplitude  required 
to  fail  an  ice  cover  of  0.5  m  thick 
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Since  the  nonlinear  term  is  important  when  the  vertical  inertia  of  water  is  significant, 
Eq.  3.42,  which  neglected  the  nonlinear  term,  is  not  strictly  correct.  Neglecting  the  friction 
and  bed  slope  and  using  the  following  dimensionless  variables 


x 


* 


x  ,*  d  ct  *  d0u 

— ,  d*  =  -,  t*  =  — ,  u  =  — 
L0  a  L0  ac 


(3.43) 


Eq.  2.9  can  be  written  in  a  nondimensional  form  as 


du*  dd*  du*  tdu* 
dF  +  dx*  +  Frdx*  +  £U  fc* 


_3V_ 

£l  dt*dx*2 


(3.44) 


where  e  =  —  and  £\  =  .  If  these  two  parameters  are  of  the  same  order,  the  nonlinear 

term  u*  should  be  considered. 


3.2.5  Group  Velocity 

The  group  velocity  plays  a  fundamental  role  in  wave  propagation  since  wave  energy  prop¬ 
agates  at  this  velocity.  In  dimensionless  form,  the  group  velocity  U  is 


ry  d/3R  „  dc 
u=-fo-c+ate 


(3.45) 


When  c  is  independent  of  a  in  the  wave  number  spectrum,  the  wave  is  not  dispersive. 
Therefore,  the  group  velocity  and  wave  celerity  are  the  same  in  kinematic,  gravity,  or 
acoustic  bands.  In  the  two  bands  where  c  is  a  function  of  a,  i.e.  dynamic  and  ice-coupled 
bands,  the  waves  are  highly  dispersive.  Thus,  the  group  velocity  and  wave  celerity  differ 
as  shown  in  Fig.  3.11  as  well  as  C.26  to  C.28.  Without  the  axial  force,  the  group  velocity 
always  exceeds  the  wave  celerity.  However,  with  a  large  enough  axial  force,  the  group 
velocity  may  approach  zero. 

When  group  velocity  varies  in  the  wave  propagation  direction  from  a  positive  value 
to  a  small  value,  there  will  be  an  accumulation  of  wave  energy.  An  infinite  wave  train 
of  initially  constant  amplitude  entering  from  a  distance,  would  attain  singular  amplitude 
where  it  encountered  a  place  of  zero  group  velocity.  From  Fig.  3.11,  The  group  velocity 
can  reach  zero  just  around  the  wave  length  2nl  over  the  entire  spectrum.  This  is  another 
explanation  for  waves  with  length  near  2ttI  can  break  ice  cover  easily. 
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Figure  3.11:  Dimensionless  wave  celerity  and  group  velocity  as  functions  of  dimensionless 
wave  number 
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Dimensionless  Wave  Celerity 


Dimensionless  Wave  Number  a 


Figure  3.12:  Enlarged  plot  of  the  wave  celerity  and  group  velocity  plot  in  Fig.  3.11 
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3.3  Summary  and  Conclusions 


In  this  chapter,  the  effect  of  an  axial  force  along  the  ice  cover  on  the  interaction  of  water 
waves  with  ice  cover  in  a  uniform  channel  is  analyzed  using  the  linear  perturbation  theory. 
The  linear  theory  is  valid  only  for  small  amplitude  waves,  i.e.  the  order  of  nonlinear  terms 
is  much  less  than  that  of  linear  terms.  The  present  analysis  shows  that  the  axial  force 
can  greatly  reduce  the  wave  celerity  in  the  ice-coupled  band  of  the  wave  number  spectrum. 
The  wave  attenuation  in  the  transition  zone  between  the  gravity  band  and  the  ice-coupled 
band  increases  with  the  axial  force.  The  fracture  of  ice  cover  will  occur  at  a  wavelength 
of  27 d  as  in  the  case  of  zero  axial  force  analyzed  by  Daly  (1995).  This  means  that  if  the 
cover  does  not  fracture  at  a  wavelength  of  27 rZ,  it  will  then  not  fracture.  This  implies  that 
the  compressibility  of  water  can  be  neglected  in  the  analysis  of  ice  cover  breakup.  The 
present  analysis  also  shows  that  the  minimum  wave  amplitude  that  is  required  to  fracture 
the  ice  cover  rapidly  approaches  zero  when  the  axial  force  approaches  the  critical  buckling 
load  of  a  floating  elastic  plate  and  the  wave  length  approaches  2irl.  The  linear  analysis  is 
therefore  valid  for  cases  with  a  large  axial  force.  The  validity  of  the  linear  theory  for  zero 
axial  force  needs  to  be  further  examined  and  will  be  discussed  in  the  Section  4.1  along  with 
the  nonlinear  analysis. 
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Chapter  4 

Nonlinear  Model 


The  linear  analysis  determined  the  asymptotic  ice  cover  response  to  propagating  waves 
by  sinusoidal  perturbations  of  various  wavelengths  and  amplitudes.  It  showed  that  waves 
with  a  wavelength  near  the  characteristic  length  of  the  ice  cover  are  most  likely  to  induce 
transverse  cracks.  Surges  that  occurred  during  the  ice  jam  release  may  make  the  nonlinear 
term  u  ^  in  Eq.2.9  significant  and  not  negligible.  No  field  data  on  water  waves  underneath 
a  river  ice  cover  is  available  for  determining  the  dominate  wave  types  which  cause  the  ice 
breakup.  To  circumvent  this  lack  of  information,  the  governing  equations  describing  the 
propagation  of  waves  under  a  floating  ice  cover  will  be  investigated  in  this  chapter  by 
considering  nonlinear  terms  of  the  water  wave  equations. 


4.1  Governing  Equations 

The  compressibility  of  water  is  not  important  when  the  ice  cover  is  separated  from  the 
banks  after  the  formation  of  the  longitudinal  cracks.  The  linear  analysis  also  showed  that 
the  compressibility  is  not  responsible  for  the  formation  of  transverse  cracks,  therefore  the 
water  is  considered  to  be  incompressible  in  the  following  analysis.  The  linear  analysis 
showed  that  friction  and  bed  slope  contribute  only  to  kinematic  waves.  Therefore,  the 
bed  and  cover  friction  can  be  assumed  to  be  small  and  the  bed  slope  is  negligible.  Thus 
the  mass  conservation  equation  2.11,  the  momentum  equation  2.8  and  the  elastic  ice  cover 
equation  2.14  can  be  written  as: 
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dd  d(ud) 

- 1-  — - — -  =  0 

dt  dx 


du  du  dd  1  dP 
«+“&+9&  +  p'&  _0 


EI  d*d  d2d  „d2d  n,  n 
1  -  u2  dx 4  +  Pin  dt2  +  N  dx2  p  ~ 


Let  u  =  u0  +  u  and  d  =  d0  +  d ,  where  it0  and  d0  are  constants,  representing  the  uniform 
flow,  then  Eqs.  4.1,  4.2  and  4.3  give 


dd  dd  du  d(ud ) 

~rrr  +  Uo~^ — I-  d0—  | - —  —  0 

dt  dx  dx  dx 


du  du  <  du  dd  1  dP'  n 

■sr+“°&  +“&  +9te  +  =0 


El  did  d2d  d2d 

T^sw+™lv+NdS-p=0  (4'6) 

Let  x  —  x  —  u0t  and  i  =  t,  thus  ^  =  ^7  and  ~  uo£r.  After  combining  Eqs. 

4.5  and  4.6,  this  set  of  equations  after  dropping  the  primes  becomes, 


dd  du  d(ud) 

m+doai  +  ~dr-° 


du  du  dd  EI  d5d  Pij].  d3d  0  d3d  2d3d .  Nd3d  n/A0. 
dt  dx  +  P  dx  (1  —  1 P)pdx5  ^  p  dt2dx  U°dtdx2  +  U°dx3  p  dx3 

It  is  convenient  to  introduce  the  following  dimensionless  parameters  to  characterize  the 

nonlinear  shallow  water  waves  based  on  a  horizontal  length  scale  L0>  a  vertical  length  scale 

d0,  and  the  wave  amplitude  o. 

c,  _  «o  _  a  x_  EI  _  /4  _  _  pi  doV  AT  _  N 

dp  (1  -v*)pSLi  Li’7  pLVN°  pgLl  (49) 
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The  meaning  of  the  first  two  parameters  are  well  known,  i.e.  Froude  number  and  dimen¬ 
sionless  wave  amplitude.  The  last  three  parameters  are  dimensionless  ice  cover  stiffness, 
ice  cover  inertia,  and  axial  force.  In  addition,  variables  in  Eqs.  4.7  and  4.8  can  also  be 
nondimensionlized  as: 


(4.10) 


in  which,  c  =  \fgda. 

In  terms  of  the  preceding  nondimensional  variables  and  parameters,  Eqs.  4.7  and  4.8 
can  be  written  in  the  nondimensional  form  as, 


dd *  du*  d(u*d*)  n 

dt*  dx*  dx* 


(4.11) 


du *  dd* 
dt*  dx* 

d3d* 

J^dt*2dx* 


tdu*  ;d5d*  ..  d3d* 

+  £U  dx*  +  Vx*  +N°d3x* 


+ 


-  2  Fr 


d3d*  2  d3  d* .  n 

dt*dx *2  dx*3’ 


(4.12) 


It  should  be  noted  that  Eqs.  4.11  and  4.12  are  similar  to  the  well-known  Boussinnesq 
equation  (Debnath  1994,  Mei  1983),  albeit  more  complicated.  In  Eq.  4.12  the  third  term 
is  the  nonlinear  term,  and  the  terms  following  it  are  dispersion  terms.  Daly’s  (1995)  linear 
analysis  showed  that  waves  with  a  wavelength  of  27r/  having  an  amplitude  in  the  order  of 
0(0.1  m)  can  cause  transverse  cracks  to  form  even  when  No  =  0.  In  such  a  case,  the  order 
of  ice  stiffness  parameter  6  =  (^)4  =  (^)4  is  0.001,  and  the  order  of  small  amplitude 
parameter  e  is  0.01  or  larger.  The  nonlinear  term  in  Eq.  4.12  is  comparable  to  the  ice 
stiffness  term,  thus  can  not  be  neglected.  The  linear  analysis  is  therefore  not  valid. 

We  look  for  a  solution,  correct  to  the  first  order  in  e  and  5,  in  the  form 


u*  =  d*  +  eQi  +  5Q2  +  0(e?  +  <J2) 


(4.13) 


where  Qi  and  Q2  are  functions  of  d*  and  its  derivatives.  Consequently,  since  d*t. 
0(e,  5),  Eqs.  4.11  and  4.12  become 

d’t.  +<-+  ecTd-,.  +£(dX.  +  ^-)  +  +  0(e\, ei)  =  0 


-d*x.  + 


(4.14) 
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(4.15) 


dQi 


dQ2  dP* 


d*r  +  d*x.  +  ed*dx .  -  £—■  +  <H— +  g~r)  +  °(£  > eS)  ~  0 

where  P*  =  +  |(g£  -  +  Fr2|^)  +  These  two  equations  must  be 

consistent  so  that  we  stipulate 


_  d*2 

<?!  -  -T 

Q2  —  ~pr 


(4-16) 


Consequently,  we  obtain  a  single  equation  for  d*  accurate  to  the  first  order,  by  neglecting 
terms  of  order  e2  and  8e  or  higher: 


d*f  +  d*x.  +  -edx.d*  +  -d 


X*X*X*X  +  X*  2"(1  +  F r)^  xm  xn 


.No. 

+  Yd 


'X*X*X* 


=  0 


(4.17) 


and  a  relation  between  u*  and  d*: 

d*2  P* 

u*  =  d*  +  e{-?j-)  +  6—  (4.18) 

It  is  noted  that  the  neglected  nonlinear  terms  in  the  ice  cover  response  equation  Eq.  2.14 
(|4)2  _  £:2(|£l)2(|4I)2  <  e2  Therefore,  the  assumption  of  £  =  is  acceptable. 

Equation  4.17  expressed  in  dimensional  form  is 

13  p  N 

“ dt  +  dx  +  -ddx  +  ~dxxxxx  +  —  d^rj(l  +  F  r)2  dxxx  +  —  dxxx  =  0  (4-19) 

c  2  d0  2  2 p  2pg 

The  first  two  terms,  dt  +  cdx,  describe  wave  evolution  at  the  shallow  water  speed  c  =  \/gda, 
the  third  term  with  coefficient  ^  represents  a  nonlinear  wave  steepening.  These  term  are 
the  same  as  those  in  the  KdV  equation  for  free-surface  shallow  water  waves.  The  rest  of  the 
terms  are  all  dispersion  terms  due  to  ice  cover  bending,  inertia  of  ice  cover,  and  the  axial 
force  along  the  cover.  Thus,  similar  to  the  KdV  equation,  Eq.  4.19  is  a  balance  between 
time  evolution,  nonlinearity  and  dispersion. 

We  now  seek  a  steady  progress  wave  solution  of  Eq.  4.19  traveling  to  the  downstream. 
The  solution  is  stationary  in  the  reference  frame  C  so  that  d  =  d( C),  C,  =  x  —  Ut  .  Substi¬ 
tuting  this  into  Eq.  4.19,  we  obtain 

(1  -  v-)i  +  £-<Ul  +  U""  +  f  4,(1  +  Frfd"  +  £-<r  =  0  (4.20) 

c  2d o  2  2  p  2pg 
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where  '  =  Integrating  with  respect  to  £,  we  get 


A  +  (1  -  -)<f  +  +  £(T  +  p-dori(l  +  Fr); V  +  2L<f  =  0  (4.21) 

v  c '  4d0  2  2p  2p0 

Where  A  is  an  integration  constant  to  be  determined. 

4.2  Solutions 

4.2.1  Solution  of  the  Simplified  Equation 

We  will  first  solve  the  following  simplified  equation  by  ignoring  the  axial  force  and  the 
inertia  of  ice  cover  terms 

A  +  (1  —  ~)d  +  +  “d""  =  0  (4.22) 

c  4  d0  2 

It  is  known  that  cnoidal  waves  with  slowly  varying  amplitude  are  observed  in  rivers  (Deb- 
nath  1994).  By  assuming  a  solution  of  the  following  form 

d(C)  =  — acn4(o:C,  k)  +  b  (4.23) 

where  a  is  the  wave  height  measured  vertically  from  trough  to  crest  and  U  is  the  wave 
celerity,  both  of  these  and  constants  a,  b  and  k  need  to  be  determined.  Since  the  net  area 
occupied  by  fluid  within  a  wavelength  above  the  mean  water  level  is  zero: 

rLw 

/  d(x,  t)dx  =  0  (4.24) 

Jo 

where  the  wavelength  as  shown  in  Appendix  B  is: 

L\v  —  2a  4K(/c)  (4.25) 

and  K (k)  is  the  complete  elliptic  integral  of  the  first  kind. 

Equation  4.22  gives 

A  +  (1  —  —)d  +  ~^rd2  =  ——a4dxxxx  (4-26) 

c  4ao  " 
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in  which  X  =  a(.  The  following  procedure  shows  that  Eq.  4.23  satisfies  Eq.  4.26.  The 
term  a4dXxxx  and  the  left-hand  terms  of  Eq.  4.26  can  be  written  as: 

a4dXXxx  =  — a[840a4A:4cn8(X)  +  1040(-2A;4  +  k2)a4cn6(X) 

+32a4(53Jk4  -  53fc2  +  8)cn4(X)  +  240a4(-2fc4  +  3A;2  -  l)cn2(A) 


+24a4(l  -  A:2)2] 

A+v-7»+ir/  -  gc„»m+[(i-f)+i|6](-„)c»<w 

a -£*+£•! 


(4.27) 


(4.28) 


The  procedure  of  deriving  (cn4)Xxxx  is  shown  in  Appendix  B.  By  equating  the  coefficients 
of  the  polynomial  of  cn  on  both  sides  of  Eq.  4.26,  the  following  five  equations  should  be 
satisfied  by  determining  values  of,  or  relationships  between,  a,  a  ,  b,  k,  A  and  U. 


I4  a4 


l4a 4  4  3  2 

— —  a840fc4  =  —  a2 
2  4d0 

al040(-2A;4  +  A;2)  =  0 


14  4  TJ  O 

a32(53A:4  —  53A:2  +  8)  =  (1 - +  - — &)(-a) 

”  2an 


2 

l4  a4 


a240(-2A:4  +  3A:2  -  1)  =  0 


l4a4 


a24(l  —  A;2)2  =  A  +  (1 


-)b+zrb2 

c  4  d0 


Equations  4.30  and  4.32  give 


k2  = 


(4.29) 

(4.30) 

(4.31) 

(4.32) 

(4.33) 

(4.34) 


Equation  4.29  yields 


and,  Eq.  4.24  yields 


in  which 


a 


140  d0l4 


,2K(*) 

/  —  acn4(y,  k)dy  +  2K(k)b  =  0 

Jo 


(4.35) 


(4.36) 


y  —  ax 


(4.37) 
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Using  k2  —  1/2  and  Eq.  B.32,  Eq.  4.36  yields 


(4.38) 


Equation  4.31  gives 


p-«‘-S5>e 


Therefore,  the  wavelength  is  given  by  Eq.  4.25  as: 


(4.39) 


(4.40) 


in  which 

K(^|)  «  1-854 

The  constant  A  is  obtained  from  Eq.  4.33  as 

23a2 

210do 


(4.41) 


(4.42) 


The  final  solution  of  Eq.  4.22  is 


d(x,  t) 


.  a  .i  (x  —  Ut) 

V140do  l 


(4.43) 


where  wave  speed  U  is  given  by  Eq.  4.39 

Equation  4.43  represents  a  train  of  periodic  waves  propagating  at  the  speed  of  U  =  (1  — 
0.1^)c  in  the  ice-covered  river  channel.  The  shape,  period,  and  wavelength  of  the  wave  are 
all  functions  of  the  wave  amplitude  a.  Equation  4.40  shows  that  the  wavelength  decreases  as 
the  wave  amplitude  increases.  Equation  4.44  shows  that  the  celerity  of  the  wave  is  smaller 
than  c.  The  effect  of  ice  cover  on  wave  speed  can  be  seen  by  comparing  Eq.  4.44  with  the 
solitary  wave  speed  U  =  c(l  +  ^),  and  cnoidal  wave  speed  U  =  c[l  +  (a  —  yj^)/2d0}  for 
free  surface  waves  (Debnath  1994).  Equation  4.43  is  an  exact  solution  of  Eq.  4.22,  which 
is  derived  under  the  condition  ^  -C  1,  for  all  positive  a. 
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Equations  4.43,  4.39  and  4.40  can  be  rewritten  as 


d 

a 


(  1  a 
^140  d0  l 


(4.44) 


U_l_J_aL 
c  10  do 


(4.45) 


and, 

^=2(i^)iK  (,/i)  (4.46) 

l  a  \  z 

A  typical  wave  profile  given  by  Eq.  4.44  is  shown  in  Fig.  4.1.  Figure  4.2  shows  profiles  of 
the  nonlinear  waves  varying  with  the  wave  amplitude  for  the  first  half  wavelength,  obtained 
from  Maple  Version  4. 


4.2.2  Application  of  Solutions  for  the  Simplified  Equation 

The  maximum  bending  stress  Sxmax  in  the  ice  cover  is  at  the  bottom  of  the  ice  sheet, 
2  =  — Eq.  3.37  becomes 

E  7}  d2d 


Sxmax  — 


1  —  V2  2  dx 2 


(4.47) 


(4.48) 


Using  Eqs.  4.43  and  B.23  with  k2  —  |,  Eq.  4.47  gives 

■W  =  -377^““2(-10cn6(X,  )/|)  +  6c"2(*,  ^/|» 

in  which,  X  =  =  (140 as  defined  in  Eq.  4.26.  The  profile  of  the  function, 

Fen  =  —  10cn6 (A,  /|)  +  6cn2(X,  /|),  is  shown  in  Fig.  4.3.  Its  minimum  value  is  -4. 
Substituting  Eq.  4.35  for  a  in  Eq.  4.4  yields 


Sxmax  — 


2 

E  r]  a* 


(4.49) 


1  -  v2 12  x/35 h 

For  an  ice  cover  with  a  maximum  flexural  strength  Sxb,  the  minimum  wave  amplitude 
which  cause  ice  cover  fracture  is 


(4.50) 
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£22*  =  ( 3  (4.52) 

do  X  E  dovj 

The  value  of  elastic  modulus  of  the  channel  ice  cover  E  ranges  from  0.4  to  9.8  GPa. 
The  value  of  Poisson’s  ratio  v  is  estimated  to  be  0.35.  The  maximum  flexural  strength  SXb 
is  assumed  to  be  0.6  MPa,  which  the  ice  cover  can  withstand  before  a  crack  forms  (Daly 
1995).  Equations  4.52,  4.39  and  4.40  are  used  to  determine  amini  and  corresponding  values 
of  U  and  Lw,  i.e.  LWmax  based  on  these  values.  Sample  values  of  the  wave  amplitude 
amin  required  to  cause  the  ice  cover  to  crack,  as  well  as  other  corresponding  parameters,  is 
shown  in  Table  4.1,  where  l  = 


Figure  4.3:  The  profile  of  the  function  — 10cn6(X,  \J\)  +  6cn2(X,  \J\) 

If  the  amplitude  of  an  incoming  wave  is  larger  than  amin  or  the  wave  length  is  less 
than  Lwmax,  the  wave  will  form  transverse  cracks  in  the  ice  cover.  The  linear  model 
gives  the  result  that  the  wave  with  wavelength  near  2x1,  i.e.  in  the  order  of  10  m  can 
generate  sufficient  amplitude  to  cause  transverse  cracks  to  form.  However,  The  nonlinear 
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Table  4.1:  Values  for  amin  and  other  corresponding  parameters 


d0  (m) 

E  (GPa) 

77  (m) 

l  (m) 

O' min  (wi) 

U  (m/s) 

L\Vmax  (ffl) 

5 

1 

0.2 

2.87 

0.475 

6.934 

65.98 

5 

1 

0.5 

5.71 

0.645 

6.910 

121.54 

5 

1 

1.0 

9.60 

0.812 

6.886 

5 

0.220 

6.969 

5 

QH| 

mm 

0.299 

6.958 

261.85 

5 

10 

1.0 

17.08 

0.377 

6.947 

415.65 

3 

1 

0.2 

2.87 

0.401 

5.350 

60.60 

3 

1 

0.5 

5.71 

0.544 

5.324 

111.62 

3 

1 

1.0 

9.60 

0.685 

5.298 

177.19 

3 

10 

0.2 

5.11 

0.186 

5.388 

130.55 

3 

10 

0.5 

10.15 

0.252 

5.376 

240.48 

3 

10 

1.0 

17.08 

0.318 

5.365 

381.73 

shown  in  Table  4.1.  This  is  more  consistent  with  the  observations  given  by  Parkinson(1982), 
Gerard  et  al.  (1984),  and  Prowse  (1986). 

The  relation  between  the  wavelength  and  wave  amplitude,  Eq.  4.40,  can  be  rearranged 


as: 

do  7 r  L\y 

This  relationship  is  shown  in  Figs.  4.4  and  4.5. 

Equations  4.23  and  B.34  give 

|-|  =  aal-cnW 
4  a5/4 


(4.52) 


-  (140do)1/4/  ^4'53^ 

Since  the  order  of  |  fj  |  is  0(0.1)  for  all  cases  in  Table  4.1,  the  assumption  used  in  deriving 
the  plate  equation,  Eq.2.14  is  justified. 
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Figure  4.4:  The  relationship  between  wavelength  and  amplitude  of  nonlinear  waves  Eq. 


Figure  4.5:  The  relationship  between  wavelength  and  amplitude  of  nonlinear  waves  in 
dimensionless  form  Eq.  4.52 


4.2.3  Solution  of  the  Complete  Equation 

Now,  we  solve  the  complete  equation  4.19.  The  idea  is  the  same  as  that  for  the  simplified 
case  discussed  in  Sec.  4.2.1.  Equation  4.21  is  rearranged  as  the  following: 

A  +  (1  -  -)d  +  i d 2  +  U"'  +  /3(f  =  0  (4.54) 

c  4cto  2 


Figure  4.6:  Characteristic  length  of  ice  cover 

where  '  =  £  =  £  =  The  characteristic  length  l  —  1S  shown  in  Fig. 

4.6.  The  parameter  j3  is  defined  as  ft  =  ^^(1  +  Ft)2  +  “2  ■  The  first  term  represents 
the  ice  inertia  term,  (3ice ,  which  is  shown  in  Fig.  4.7.  The  second  term  is  the  axial  force 
Paxiai ■  If  N  equals  the  maximum  axial  force  defined  in  Eq.  3.33,  then  fdaxiai  —  1- 
Assuming  the  solution  is  in  the  following  form 

d(  C)  =  — acn4(a:£,  k )  +  6cn2(a£,  k)  +  g  (4.55) 

i.e. 

./  .  i,  X—Ut  ..  ,  9/  x  —  Ut 

d(x,  t )  =  — acn  (a — - — ,  k )  +  6cn  (a — - — ,  k)  +  g 

l  V 
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where  a,  a,  b,  U,  k  and  g  are  to  be  determined.  The  coefficient  a,  b,  and  a  are  different 
from  those  in  Eq.  4.23.  The  wave  height  in  Eq.4.55  is  a  combination  of  a  and  b. 

Equation  4.54  can  be  rewritten  as  the  following  by  letting  X  =  a£. 

-[A  +  (1  -  -)d  +  ^-d2]  =  \a4dXxxx  +  I3a2dxx  (4-56) 

Equation  4.56  is  similar  to  Eq.4.26,  except  for  the  last  term.  The  right-hand  side  terms  in 
Eq.4.56  can  be  obtained  as  the  following  by  substituting  Eq.4.55  into  Eq.4.54.  Expressions 
for  cn?xxxx  and  ciixx  are  derived  in  Appendix  B. 

-a4dxxxx  —  —  a840a4-fc4cn8(A) 

+[— al040(— 2A;4  +  k2)a 4\  +  68  x  ^W^jcn6^) 

z  z 

+[— a32(53A;4  -  53A;2  +  8)^a4  +  68  x  15(— 2A:4  +  k2)^a4} cn4(X) 

+[-a24oJa4(-2A;4  +  3 k2  -  1)  +  6^48(17A:4  -  17fc2  +  2)]cn2(AT) 
z  z 

+ (-a)24^a4(l  -  k2)2  +  bla48{-2k4  +  3 k2  -  1)  (4.57) 

z  z 

/3a2dxx  =  a2[/?(20afc2)cn6(A) 

+  (-ap  16(— 1  +  2fc2)  4-  bp{-6k2))  cn4(X) 

+  {-apl2(l  -  k2)  +  6/?4(— 1  +  2 A:2))  cn2(X) 

+bp2{\  -  A;2)]  (4.58) 

The  left-hand  side  terms  can  be  obtained  as: 

A  +  (\--)4+%  =  ^cn!(X)  +  4-(-o)<«i>8W  +  [(l  +  2§-(62  +  2(-o)9)lcn4m 

c  4  d0  4  d0  2  d0  c  4d0 

+[(l--)i>  +  +  (A  +  (  !-%+¥■)  (4.59) 

c  4ao  c 

By  equating  the  coefficients  of  polynomials  of  cn  on  the  two  sides  of  Eq.  4.56,  The 
following  five  conditions  are  obtained  for  determining  a,  a,  6,  U,  k  and  g,  among  which 
only  one  is  an  independent  parameter. 

— ^-a2  =  —  a840a4^A;4 
4d0  2 
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~2  (~a)b  =  (-a)m0{-2k4  +  k2)aAl 
4do  * 

+a4\b8  x  15A:4  +  (20 apk2)a2 

z 


(4.60) 


.[(1  _  — )(— a)  +  -^-(62  +  2(-a)5)]  =  (-a)32a4i(53A:4  -  53A;2  +  8)  +  b\a48  x  15(-2A:4  +  k 2) 
c  4ao  2  ^ 

+(— a)/?16(— 1  +  2k2)a2  +  bp(-6k2)a2  (4.61) 

8(17/b4  -  17& 

(4.62) 


_[(1_^)6  + JL2fy)]  =  (— a)240a4^(— 2A:4  +  3A:2  —  1)  +  i>0!4^8(17A:4  —  17A:2  +  2) 
c  4do  2  2 

+(-a)/?12(l  -  A:2) a2  +  b(34(- 1  +  2A;2)a2 


+-j-/?2(l  -  A;2)a2 
do 

The  integration  constant  A  needs  to  be  determined.  Since  the  net  area  occupied  by 
fluid  within  a  wavelength  above  the  mean  water  depth  do  is  zero: 

rLw 


b  A 


(4.63) 


rLw 

/  d(x,  t)dx  =  0 
Jo 


(4.64) 


which  implies 


J^W[—acn4  ^a—  k^j  +  6cn2  ^*(x  kj  +  g]dx  =  0  (4.65) 


where  the  wavelength  (see  Appendix  B) 


2/ 

Lw  =  — K(A;) 

a 


Let 


y  =  a 


(x  -  Ut) 

l 


Equation  4.65  yields 

,2K(*) 


rl  K(fc)  ,  ,2K(fc)  „ 

—a  /  cn 4(y,  k)dy  +  b  cn 2(y,  k)dy  +  2K(A;)p  =  0 

Jo  Jo 

Using  Eqs.  B.31  and  B.32,  Eq.  4.68  yields 

9  _  f2(2fc2  -  1)  b\  ( 

a  \  3k2  a)  \  K(k)k 2 


g  _  (2{2k2  -  1)  b\(  E (fc)  l-k2\  1  —  A:2 


k2 


3k2 


(4.66) 


(4.67) 


(4.68) 


(4.69) 
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The  six  equations,  4.59  to  4.63,  and  4.69,  are  used  to  determine  the  constant  A  and  five  of 

the  six  parameters  a,  a,  b  ,U,  k  and  g.  Now,  the  parameter  k,  which  varies  between  0  and 

1,  is  selected  as  the  independent  variable.  Equations  4.59  to  4.62  give 

j  =  560 k4a4  (4.70) 

i  =  <4-71> 

-[(1  -  +  |f-]  =  (l6(53fc4  -  S3 k2  +  8)  +  (40  -  ^^)(1  -  2fc2)2)  a* 

+<-2»  +  +  W*  ~  2  s  0  t  f  (4  72) 

+ »*» - 1) + cm  -  * V)  2{_J+?)a,  _  a 

+4(17k4  -  17fc2  +  2 )a4  +  4(-l  +  2 k2)/3a2  (4.73) 

After  combining  Eqs.  4.72  and  4.73  ,  we  obtain  the  following  equation  for  j- 

C3(j)3  +  C2(j)2  +  Ct(j)  +  Co  =  0  (4.74) 

where 


C3  =  3  x  120( — 2A;4  +  3A:2  —  l)k2  +  8(17&4  —  17fc2  4-  2){—2k2  +  1) 

— 2(— 2 k2  +  1)  (l6(53fc4  -  53fc2  +  8)  +  (40  -  ^yi*)(l  -  2fc2)2) 

C2  =  [3fc212(l  -  A:2)  —  ^-(llk,4  -  YJk2  +  2)  -  8(-l  +  2fc2)2 

lo 


'2(_20+3lfl3)(“2‘:2  +  1)! 


(l6(53 k4  -  53k?  +  8)  +  (40  -  7°  *  8)(1  -  2 fc2)2)  /13] 


„  r 4  /  on  500  ,  1  4  x  31 

Ci  -  l13  +  (  +  3  x  13^  13  +  3  x  132 


](-2  k2  + 1) 


Co  =  - 


2  x  31 

3  x  13s 


Thus  is  a  function  of  k  which  can  be  obtained  from  Eq.  4.74.  Furthermore,  should 
be  the  positive  real  roots  of  Eq.  4.74.  Equations  4.70  and  4.71  yield 


b  _  1  _  2(1  -  2 fc2) 

a  ~  39k2  f  3  k2 
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Thus  Eq.  4.69  yields 


g  -if  E  (k)  1  —  k2\  1  —  k2 

a~39Pf{  K(k)k2  k2  )  +  3k2 


(4.75) 


This  shows  that  *  is  a  function  of  k  and  is  independent  on  /?.  Eq.  4.72  is  rearranged  as 
the  following, 


3  a  g 
2  do  a 


+  Q 


(4.76) 


where  Q  is  the  right-hand  side  term  of  Eq.  4.72  and  is  a  function  of  k  and  /?.  Since  J  can 
be  replaced  by  Eq.  4.75,  U  is  a  function  of  k  and  /?. 

The  following  procedure  is  used  to  establish  curves  for  a,  b,  U,  g  and  a  as  functions  of 
k  for  a  selected  value  of  the  physical  parameter  /?: 

1.  Select  a  physical  parameter  /?; 

2.  Choose  a  value  of  k,  which  is  in  the  range  of  0  to  1; 

3.  Calculate  y  and  a:  from  Eq.  4.74; 

4.  Calculate  coefficient  a  from  Eq.  4.70; 

5.  Calculate  b ,  g  and  U  from  Eqs.  4.71,  4.75,  and  4.76,  respectively,  and  the  wave 
profile  d  from  Eq.  4.55; 

6.  The  wave  length  Lw  can  be  calculated  using  Eq.  B.9,  i.e. 


(4.77) 


Figures  4.8  to  4.12  show  the  dimensionless  parameters  a,  b,  a,  g  and  U  as  functions  of  k. 
Similar  plots  for  a ,  a,  b  and  g  in  dimension  are  shown  in  Appendix  D.  These  figures  show 
that,  when  k  ->■  0,  thus  a  -f  a  constant  which  depends  on  the  flow  and  ice  cover  condition 
l  and  /?,  and  a  -4  0,  b  -»  0.  The  result  approaches  that  of  linear  equation.  When  k  ->  yj\, 
a  ,  a  and  b  all  approach  positive  infinity.  In  this  case,  the  result  approaches  that  of  the 
simplified  equation. 
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Table  4.2:  Solutions  of  the  complete  equation  for  Lw  —  381.7  m 


p 

a(m) 

b(m ) 

k 

a 

9 

0.0079 

0.323 

0.01134 

0.7126 

0.16524 

0.3164 

-0.01146 

0.0130 

0.340 

0.0192 

0.7160 

0.16660 

0.3062 

-0.01245 

0.3130 

0.330 

0.2362 

0.8382 

0.18593 

0.00047 

-0.0678 

Given  a  wavelength  Lw  —  381.7  m,  and  let  others  parameters  Fr  =  0.3,  d0  =  3  m, 
E  -  10  GPa,  and  77  =  1  m  be  the  same  as  those  of  the  last  case  in  the  Table  4.1  for  the 
simplied  solution,  solutions  of  Eq.  4.54,  which  obtained  from  Figs.  4.8  to  4.12,  are  shown 
in  Table  4.2.  Figure  4.13  shows  the  profiles  of  these  nonlinear  waves.  When  /?  =  0,  this 
returns  to  the  simple  case.  When  only  ice  inertia  is  considered,  /3  =  0.0079,  i.e.  Case  1  in 
Table  4.2.  If  there  is  an  axial  force  in  ice  cover,  say  that  is  Case  3  in  Table  4.2. 


Figure  4.12:  Variation  of  parameter  U  with  k 
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Figure  4.13:  Profiles  of  nonlinear  waves 


4.2.4  Application 

The  maximum  bending  stress  in  the  ice  cover  Sxmax  is  at  the  bottom  of  the  ice  sheet, 


—  _2 
~  2' 


,  _  E  r]d2d 

>Xmax  -  l_v22dx2 


(4.78) 


Based  on  Eqs.  4.55,  B.23  and  B.20,  can  be  obtained  as  Eq.4.79,  after  setting  t  =  0, 


^-[20A;2acn6(a:y,  k)  +  (16(— 1  +  2 k2)(-a)  —  6k2b)cni(aj,  k ) 

+[12(1  -  k2)(—a)  +  2(— 2  +  4&2)6]cn2(ay,  k)  +  26(1  -  k2)}  (4.79) 

Therefore,  Eq.4.78  gives 

Sxmax  =  y^|^-[20fc2acn6(ay,  k)  +  (16(— 1  +  2 k2)(-a)  -  6fc26)cn4(o!y,  k) 

+[12(1  -  k2)(—a)  +  2(— 2  +  4fc2)6]cn2(ay,  k)  +  26(1  -  k2)]  (4.80) 


dx2 
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y 


Figure  4.14:  Variation  of  cnoidal  functions  with  modulus  k  =  \ 

The  profiles  of  the  functions,  cn2(x,  \J\),  cn4(x,  yj\),  and  cn6(x,  are  shown  in  Fig. 
4.14.  Figure  4.14  shows  that  their  wave  lengths  are  the  same,  as  stated  by  Eq.  B.9.  Let 

G{x,Lw)  =  ^■[20k2acn6(aj,k)  +  (16(-l  +  2k2){-a) -6k2b)cn4(aj,k) 

+[12(1  -  k2)(-a)  +  2(— 2  +  4fc2)6]cn2(<*y,  k )  +  26(1  -  A:2)]  (4.81) 

Equation  4.81  gives 

S*moI  =  7-^?C?m0I(Lw)  (4-82) 

The  maximum  value  Gmax{Lw)  is  difficult  to  obtain  explicitly.  It  can  be  obtained  nu¬ 
merically  by  varying  the  value  of  x  over  one  wave  length.  If  the  value  of  Sxmax  is  found 
to  be  larger  than  the  maximum  flexural  strength  Sxb,  then  the  ice  cover  will  crack.  The 
longest  wave,  Lwmax >  to  cause  cracking  can  be  then  found.  In  summary,  the  procedure  for 
determining  LWmax  is: 

1.  Select  a  k  (  This  fixed  Lw  ); 
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Table  4.3:  Solutions  of  complete  equation  for  breakuping  ice  cover 


p 

Lw(m) 

a(m) 

b(m) 

k 

a 

9 

381.7 

0.3180 

0.0 

0.3333 

-0.0106 

0.013 

387.5 

0.3230 

0.0186 

0.7162 

0.1645 

0.3056 

-0.0120 

0.113 

414.9 

0.4310 

0.1836 

0.1617 

-0.0237 

0.213 

440.8 

0.3598 

IWBHB 

-0.0349 

0.313 

442.4 

imu 

0.4793 

0.8432 

-0.0478 

-0.0504 

2.  Vary  x  to  find  the  maximum  value  of  Sxmax  —  Lw)  for  this  k  from 

Eq.  4.82; 

3.  See  if  this  matches  or  exceeds  Sxb ■  If  not,  try  another  k,  and  repeat  steps  1  to  3. 
Using  the  parameters  of  the  last  case  in  the  Table  4.1,  with  Fr  =  0.3,  waves  which  cause 

the  ice  cover  to  fracture  are  shown  in  Table  4.3.  The  profiles  of  these  waves  are  shown  in 
Fig.  4.15. 

If  an  axial  force  is  absent,  the  effect  of  inertia  of  ice  cover  is  insignificant  as  shown 
in  Table  4.2  and  Fig.  4.13.  If  the  axial  force  is  large  enough,  it  can  change  the  shape 
of  the  wave.  Furthermore,  the  wave  height  required  to  cause  the  ice  cover  to  fracture  is 
significantly  reduced,  as  shown  in  Fig.4.15. 

Eqs.  4.55,  B.33  and  B.34  give 


=  T  l-4“4W  +  *s-’WI 

<  y max{a,  6)(|  ^cn4(x)  |  +  |  ^cn2(x)  |) 

<  ^^—-max(a,  b )  (4  +  2) 

Lw 


< 


12K  (k) 
Lw 


max(a,  b) 


Since  a  <  1  m,  b  <  1  m,  K (k)  <  3  as  k  <  9  and  the  order  of  wavelength  is  0(100)  m,  and 
the  order  of  |  |  is  0(0.1),  the  assumption  used  in  deriving  the  plate  equation  Eq.  2.14  is 

again  justified. 
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d(m) 


6 


4.3 


Summary 


In  this  chapter,  a  nonlinear  analysis  is  carried  out  for  the  interaction  of  shallow  water 
waves  with  a  floating  ice  cover  in  a  uniform  channel.  By  retaining  the  nonlinear  terms 
of  order  e  and  6  (dimensionless  wave  amplitude  and  dimensionless  ice  cover  stiffness),  a 
nonlinear  wave  equation  similar  to  the  form  of  the  Korteweg-de  Vries  (KdV)  equation  is 
obtained.  The  nonlinear  periodic  wave  solutions  similar  to  the  cnoidal  wave  are  derived 
for  this  nonlinear  equation.  The  solution  for  a  simplified  form  of  this  equation,  neglecting 
the  axial  force  and  the  ice  inertia  term,  shows  that  the  celerity  of  the  shallow  water  wave 
is  slightly  reduced  by  the  existence  of  the  ice  cover.  The  minimum  wave  height  that  is 
required  to  fracture  the  ice  cover  is  typically  in  the  range  of  0.2  to  0.8  m,  depending  on  the 
ice  cover  thickness  and  strength.  The  corresponding  wavelength,  i.e.  the  distance  between 
the  transverse  cracks  that  formed  by  the  propagating  water  waves,  varies  from  50  to  400  m. 
For  typical  ice  cover  conditions  during  the  spring  breakup  period,  the  range  is  in  the  order 
of  50  to  150  m.  These  results  agree  closely  with  the  existing  field  observations  by  Parkinson 
(1982),  Gerard  (1984),  and  Prowse  (1986).  The  solution  of  the  complete  equation  shows 
that  the  ice  inertia  term  has  an  insignificant  effect.  However,  an  axial  force  could  change 
the  wave  profile  and  reduce  the  wave  height  if  the  wavelength  is  fixed.  Furthermore,  the 
axial  force  reduces  the  minimum  wave  height  that  is  required  to  cause  ice  cover  fracture. 
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Chapter  5 


Conclusions 


In  this  study  the  interaction  of  shallow  water  waves  with  an  ice  cover  in  a  river  channel 
is  analyzed  using  both  linear  and  nonlinear  theories.  The  linear  analysis  shows  that  the 
linear  theory  developed  by  Daly(1995),  which  neglected  the  effect  of  the  axial  force  acting 
along  the  ice  cover,  is  limited  to  small  amplitude  waves  and  can  not  be  used  to  determined 
the  condition  for  ice  cover  fracture.  The  inclusion  of  the  axial  force  reduced  the  magnitude 
of  the  wave  amplitude  that  is  required  to  fracture  an  ice  cover  so  that  validity  of  the  linear 
analysis  is  satisfied.  However,  the  results  obtained  from  the  linear  analysis  do  not  agree 
with  the  existing  field  observations.  A  further  analysis  by  including  the  nonlinear  terms 
in  the  water  wave  equation  is  needed.  In  the  nonlinear  analysis,  a  nonlinear  equation  is 
obtained  by  retaining  the  first  order  terms  in  the  nonlinear  governing  equations  of  the 
water  wave  and  floating  ice  cover.  The  periodic  wave  solutions  similar  to  cnoidal  waves 
are  obtained.  These  solutions  show  that  the  celerity  of  the  water  wave  is  slightly  reduced 
by  the  ice  cover.  The  minimum  wave  amplitude  that  is  required  to  fracture  the  ice  cover 
during  the  spring  breakup  period  is  in  the  order  of  0.5  m.  The  corresponding  wavelength, 
i.e.  the  spacing  between  transverse  cracks  induced  by  the  water  wave,  is  in  the  order  of  50 
to  150  m.  These  results  agree  closely  with  the  existing  field  observations  on  the  formation 
and  propagation  of  breaking  front  in  the  field.  This  study  provided  a  plausible  explanation 
of  the  formation  of  the  transverse  cracks  during  the  passage  of  a  river  wave,  which  is  the  key 
mechanism  of  the  initiation  of  breakup  ice  runs  and  ice  jams.  Further  field  and  laboratory 
studies  should  be  carried  out  to  validate  the  present  finding. 
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Appendix  A 
Notations 


The  following  symbols  are  used  in  this  report: 

A  =  cross-sectional  area  of  channel; 

c  =  wave  celerity; 

cn(x,  k)  =a  Jocobian  function; 

d  =  depth  of  channel  from  bottom  to  ice  cover; 

do  =  uniform  flow  depth; 

d  —  small  perturbation  to  depth; 

E  =  elastic  modulus  of  ice; 

Ew  =  fluid  bulk  modulus  of  elasticity; 

Fr  =  channel  Froude  number; 
g  =  acceleration  of  gravity; 

H  =  piezomotric  head; 

H0  =  piezomotric  head  under  uniform  flow  condition; 
H’  =  perturbation  to  piezomotric  head  Ho] 

I  =  moment  of  inertia  of  ice  cover; 

L  =  wavelength; 

Lo  =  longitudinal  length  scale; 

Lw  =  nonlinear  wave  length  ; 
l  =  characteristic  length  of  ice  cover; 

^o  =  x  » 
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k  =  modulus  of  Jacobian  functions; 

K(k )  =  complete  elliptic  integral  of  the  first  kind; 

N  =  an  axial  force  in  ice  cover; 

N0  =  a  nondimensional  axial  force; 

P9^0 

P'  =  pressure  deviation  from  hydrostatic  beneath  ice  cover; 
Sf  =  friction  slope; 

So  =  channel  bottom  slope; 

Sx  =  stress  in  ice  cover; 

T  =  wave  period; 
t  =  time; 

u  =  mean  flow  velocity; 

U  =  celerity  of  a  nonlinear  wave; 

U  =  group  velocity; 

u0  =  mean  flow  velocity  under  steady  flow; 
u  =small  perturbation  to  flow  velocity; 
x  =  longitudinal  distance; 
z  =  elevation  of  channel  above  datum; 
a  —  nondimensional  fluid  bulk  modulus; 

P  =  nondimensional  complex  wave  propagation  factor; 

Pi  =  imaginary  part  of  propagation  factor; 

Pn  =  real  part  of  propagation  factor; 

P*  =  rescaled  propagation  factor; 

5  =  natural  decrement;  ice  characteristic  of  ice  stiffness 

rj  =  ice  cover  thickness; 

r)o  =  nondimensional  ice  cover  thickness; 

7  =  nondimensional  ice  inertial; 


v  =  Poisson’s  ratio  of  ice  cover; 
p  =  fluid  density; 
p  =  ice  density; 
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a  =  nondimensional  wave  number; 
r  =  average  shear  stress  acting  on  bed  and  ice  cover; 
t0  =  average  shear  stress  under  steady  flow  conditions;  and 
t  =  small  perturbation  to  the  average  shear  stress, 
e  =  small  amplitude  parameter. 


Appendix  B 

Elliptical  Functions 


The  detail  of  elliptical  functions  is  available  in  books  by  Lawden  (1989),  Byrd  and  Friedman 
(1971)  and  others.  A  brief  introduction  to  elliptic  functions  and  theirs  integrals  is  presented 
in  this  Appendix. 

Consider  the  integral 

X(y)  =  [V  -  =  (B.l) 

^  1  —  A:2sin2(t) 

where  the  parameter  k,  is  taken  such  that  0  <  k  <  1.  Eq.  B.l  may  be  compared  to 

•M-jfTCT  (B'2) 

When  we  use  t  =  sin(0)  so  that  w  =  arcsin(u)  or  v  =  sin(tu).  This  led  Jacobi  to  define  a 
new  pair  of  inverse  functions  from  Eq.  B.l  as  elliptic  functions: 


sn(x,  k)  =  sin(y),  cn(a:,  k)  =  cos(y)  (B.3) 

where  k  is  the  modulus.  In  addition,  a  third  elliptic  function  is  defined  as: 


dn(x,  k )  =  \J\  —  A;2sin2(y) 
The  period  of  cn(x,  A:)  and  sn(x,  k)  can  be  written  as 


(B.4) 
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(B.5) 


d6 

(1  —  k2  sin2#)  2 


4/f _ * _ _ 

Jo  (1  —  k2sm29)i 


Figure  B.l:  The  elliptic  integral  function  of  the  first  kind 


The  latter  equation  gives  the  complete  elliptic  integral  of  the  first  kind 


d6 

(1  —  A;2sin20)2 


(B.6) 


The  values  of  K(k)  in  Fig.  B.l  are  available  in  a  mathematical  handbook,  (Abramowitz 
and  Stegun  1972). 

The  two  special  cases  k  =  0,  and  1  enable  integrals  (B.l)  and  functions  (B.3)  to  be 
reduced  to  elementary  functions.  If  k  =  0  then, 


71 


x  =  y, 


cn(x,  0)  =  cos(y)  =  cos(x) 


(B.7) 


It  is  immediately  clear  that  K(0)  =  |.  If  k  —  1  the  integral  can  be  evaluated  to  yield 


x  =  arcsech(cos(y)),  cn(x,  1)  =  sech(x) 


(B.8) 


Both  cn(x,k)  and  sn(x,k)  are  periodic  functions  for  0  <  k  <  1. 

The  wavelength  of  cn(ax,  k)  is  4a-1K(A;).  The  wavelength  of  cnoidal  waves  cn2(ax,  k), 
cn4(ax,  k )  and  cn6(ax,  k)  are 


2o-1K(&) 


(B.9) 


Using  cn(x)  as  the  abbreviation  of  cn(x,  k),  the  following  equations  can  easily  be  ob¬ 
tained  from  B.3  and  B.4 


Assuming 


then 


sn2(x)  +  cn2(x)  =  1 

(B.10) 

dn2(x)  +  fc2sn2(x)  =  1 

(B.ll) 

1 

dn2(x)  —  &2cn2(x)  =  1  —  k2 

(B.12) 

.  S“W 

=  -p— cn  =  v/l  —  fc2sin2(y)  — cos(y) 
ax  ay  v  ay 

(B.13) 

=  — sn(x)dn(x) 

(B.14) 

=  cn(x)dn(x) 

(B.15) 

=  — A;2sn(x)cn(x) 

(B.16) 

y(x)  =  cn2(x) 

(B.17) 

y  (x)  =  — 2cn(x)sn(x)dn(x) 

(B.18) 
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y'2{x)  =  4  ((1  -  A:2)cn2(x)  +  (-1  +  2A;2)cn4(:r)  -  fc2cn6) 
cn2(x )  =  2  (— cn2(x)dn2(a;)  +  sn2(a;)dn2(a;)  +  fc2sn2(x)cn2(a;)) 

=  2  (-3A;2cn4(:r)  +  (—2  +  4A;2)cn2(x)  +  1  -  A;2) 
y"{x)  =  8[cn(z)sn(a;)dn3(x)  +  k2  sn(x)dn(rc)cn3(x) 

— fc2cn  (a:)  dn  (x)  sn3  (a:)] 

-j-jCti2(x)  =  8[15fc4cn6(a;)  +  (— 30/c4  +  15/:2)cn4(x) 
ax 

+(17fc4  -  17&2  +  2)cn2(z)  +  {-2k4  +  3 k2  -  1)] 

d2  4  /  \  d2  2 
^cn(x)  =  —2y 

=  2  y'2  +  2  yy" 


t  /!/  _  //// 


-cn4(:r)  =  6  y"  +8  yy  +  2yy 


840&4cn8(x)  +  1040(— 2A:4  +  fc2)cn6(x) 
+32(53A:4  -  53&2  +  8)cn4(x)  +  240(-2fc4 
+Sk2  -  l)cn2(rc)  +  24(1  -  k2)2 


(B-19) 


(B.20) 


(B.21) 


(B.22) 


=  —  20fc2cn6(x)  +  16(-1  +  2k2)cn4(x)  +  12(1  —  k2)cn2(x)  (B.23) 


(B.24) 


The  complete  integral  of  the  second  kind  is  denoted  by  E  and  is  defined  by 

/-K(fc)  /-f  i 


E  (k)  =  [  dn 2udu  =  [*  (1  —  &2sin20)2< 

Jo  Jo 


(B.25) 


The  values  of  E(/c)  are  available  in  (Abramowitz  and  Stegun  1972).  Some  special  relations 
and  values  are  listed  as  following: 


sn(0,  A;)  =  0 


sn(u  +  2K(fc),  k) 


-sn(u,  k) 


cn(0,  k)  =  1 

cn(«  +  2K(k),  k)  =  — cn (u,k) 
,2K(fc)  „ 

/  dn  2{x,k)dx  =  2E(fc) 

Jo 


(B.26) 

(B.27) 

(B.28) 

(B.29) 

(B.30) 


!(x  ,k)dx  =  jf  (\-Ldn  2(x,k)-^—^-)dx 
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Eq.  B.18  yields 


2E(fc)  1  -  k2 


2K(k) 


k 2  k2 

/o  '  3*!  3fc2 


1  d2  2 


cn2(x,  A;))efa; 


(B.31) 


6  k2  dx2 

2(2k2  —  1)  y2K(*)  1-fc2 

— sWo  cn(x,t)<fa  +  -sr2K(*) 

-^j(-2cn(x,A;)sn(i,fc)  dn(a:,fc))  |2Kw 

cn?(x,  k)ix  +  ^2KW  (B.32) 


den2  (a:) 
dx 


l<2 


(B.33) 


dcn4(x) 

dx 


|=|  2cn2(x) 


dcn2(x) 

dx 


l<4 


(B.34) 


To  order  to  check  Eq.  B.24,  The  two  special  cases  i.e  fc  =  0  and  A;  =  1  are  considered. 


— cn4(x,0)  =  256cn4(x,0)  +  (— 240)cn2(x,0)  +  24 


(B.35) 


^jcn4(x,  1)  =  840cn8(x,  1)  +  (-1040)cn6(x,  1)  +  256cn4(x,  1) 

Also,  we  can  obtain  from  the  MAPLE  4: 
d4 


(B.36) 


dx4 


cos4(x)  =  24sin4(x)  -  192cos2(x)  +  40cos4(x) 


=  256cos4(x)  +  (-240)cos2(x)  +  24  (B.37) 


sech4(x)  =  256sech4(x)tanh4(x)  —  528sech4(x)tanh2(x)(l  -  tanh2(x)) 
dx4 

+56  sech4(x)(l  —  tanh2(x))2 
=  256sech4(x)(l  —  sech2(x))2 

— 528sech4(x)sech6(x)(l  —  sech2(x))  +  56sech8(x) 

=  840sech8(x)  +  (— 1040)sech6(x)  +  256sech4(x)  (B.38) 


74 


Therefore, 


d4 


dx4 

d*__ 

dx4 


sech4(x) 


cos4(x) 


d?_ 
dx 4 


cn4(x, 1) 


d4 

dx4 


cn4(x,0) 


(B.39) 

(B.40) 
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Appendix  C 

Figures  in  Linearized  Models 


Figure  C.l:  Effect  of  elastic  modulus  on  the  influence  of  axial  force  on  wave  celerity 
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Dimensionless  Wave  Celerity  Q  Dimensionless  Wave  Celerity 


2:  Effect  of  elastic  modulus  on  the  influence  of  axial  force  on  wave  celerity 


Figure  C.3:  Effect  of  the  magnitude  of  the  axial  force  on  wave  celerity 
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Figure  C.8:  Effect  of  the  magnitude  of  the  axial  force  on  wave  celerity 


Figure  C.9:  Effect  of  the  magnitude  of  the  axial  force  on  wave  celerity 
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Figure  C.10:  Effect  of  the  magnitude  of  the  axial  force  on  wave  celerity 


Figure  C.ll:  Effect  of  the  magnitude  of  the  axial  force  on  wave  celerity 
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Decrement  (-8)  Decrement  (-8) 
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Nondimensional  Wave  Number  a 


Figure  C.12:  Effect  of  axial  forces  on  natural  decrement 


Nondimensional  Wave  Number  o 


Figure  C.13:  Effect  of  axial  forces  on  natural  decrement 


Figure  C.14:  Dimensionless  free- wave  celerity  Ch  and  dimensionless  wave  celerity  c  in 
ice-covered  channel  as  functions  of  dimensionless  wave  number 
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Dimensionless  Wave  Celerity  Dimensionless  Wave  Celerity 


Figure  C.15:  Dimensionless  free-wave  celerity  Ch  and  dimensionless  wave  celerity  c  in 
ice-covered  channel  as  functions  of  dimensionless  wave  number 
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Figure  C.16:  Dimensionless  free-wave  celerity  Cu  and  dimensionless  wave  celerity  c  in 
ice-covered  channel  as  functions  of  dimensionless  wave  number 


Figure  C.17:  Dimensionless  free- wave  celerity  Ch  and  dimensionless  wave  celerity  c  in 
ice-covered  channel  as  functions  of  dimensionless  wave  number 
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Figure  C.18:  Dimensionless  free-wave  celerity  Ch  and  dimensionless  wave  celerity  c  in 
ice-covered  channel  as  functions  of  dimensionless  wave  number 
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Minmum  wave  amplitude  (m)  Q  Minmum  wave  amplitude  (m) 


.19:  Minimum  wave  amplitude  (m)  required  to  fail  an  ice  cover  of  0.5  thick 
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Minmum  wave  amplitude  (m)  O  Minmum  wave  amplitude  (m) 


Minmum  wave  amplitude  (m) 
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Dimensionless  Wave  Celerity  Dimensionless  Wave  Celerity 


Dimensionless  Wave  Celerity  Dimensionless  Wave  Celerity 
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Figure  C.25:  (a)Effect  of  Boussinesq’s  term:  without  axial  force;  (b)  with  an  axial  force. 
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Figure  C.26:  Dimensionless  wave  celerity  and  group  velocity  as  functions  of  dimensionless 
wave  number 
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Figure  C.27:  Dimensionless  wave  celerity  and  group  velocity  as  functions  of  dimensionless 
wave  number 
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Dimensionless  Wave  Celerity  Dimensionless  Wave  Celerity 


Figure  C.28:  Dimensionless  wave  celerity  and  group  velocity  as  functions  of  dimensionless 
wave  number 


Appendix  D 

Figures  in  Nonlinear  Model 


Figure  D.l: 
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Figure  D.12: 


